!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth 
! Produced at the School of Marine Science & Technology 
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For 
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
! 
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu 
! The full copyright notice is contained in the file COPYRIGHT located in the 
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
! AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.  
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/

MODULE MOD_NORTHPOLE
   USE ALL_VARS 
   USE MOD_SPHERICAL
   USE MOD_PAR
   USE MOD_UTILS

   IMPLICIT NONE
   SAVE
   INTEGER :: NODE_NORTHPOLE            !Node index at the north pole point
   INTEGER :: MP,NP,NPE,NPCV
   INTEGER, ALLOCATABLE :: NODE_NORTHAREA(:)
   INTEGER, ALLOCATABLE :: CELL_NORTHAREA(:)   
   INTEGER, ALLOCATABLE :: NPEDGE_LST(:)   
   INTEGER, ALLOCATABLE :: NCEDGE_LST(:)   
   INTEGER, ALLOCATABLE :: MP_LST(:),NP_LST(:)   
   REAL(DP), ALLOCATABLE :: A1U_XY(:,:),A2U_XY(:,:)
   REAL(DP), ALLOCATABLE :: AW0_XY(:,:),AWX_XY(:,:),AWY_XY(:,:)

#  if defined (SEMI_IMPLICIT)
   REAL(SP), ALLOCATABLE :: XFLUX3_NP(:,:), YFLUX3_NP(:,:)
   REAL(SP), ALLOCATABLE :: UBETA_NP(:,:),VBETA_NP(:,:)
#  endif

   CONTAINS
!==============================================================================|
     SUBROUTINE FIND_NORTHPOLE
     IMPLICIT NONE
     INTEGER :: I,ITMP,NODE_NORTHPOLE_GL,IERR,NDOM
     INTEGER,ALLOCATABLE :: TMP(:)
     
     NODE_NORTHPOLE = 0
!     NODE_NORTHPOLE_GL = 0
     NP = 0
     MP = 0
     NPE = 0
     NPCV = 0

     NDOM=0
     DO I = 1,MT
!        IF(ABS(VY(I)-90.0_SP) .LE. (NEAREST(90.0_SP,1.0_sp) - 90.0_SP))THEN
        IF(ABS(VY(I)-90.0_SP) .LE. 10E-4)THEN
           NODE_NORTHPOLE = I
           ndom=ndom+1
        END IF
     END DO
     
     IF (nDOM .GT. 1) CALL Fatal_Error('Found more than one north pole node in the grid',&
          &'This should never happen, TGE should crash first?')


     ! SINCE SPHERICAL NOW ALWAYS LOOKS FOR A NORTH POLE, DO NOT CALL
     ! FATAL_ERROR IF WE DON'T FIND A NORTH POLE!

     NODE_NORTHPOLE_GL = NGID_X(NODE_NORTHPOLE)


# if defined(MULTIPROCESSOR)
     IF (PAR) THEN
        ! GATHER THE NORTH POLE NODE GLOBAL ID TO ALL PROCS AND COMPARE RESULTS
        ALLOCATE(TMP(NPROCS)); tmp = 0
        CALL MPI_ALLGATHER(NODE_NORTHPOLE_GL,1,MPI_INTEGER,TMP,1,MPI_INTEGER,MPI_FVCOM_GROUP,IERR)

        NDOM=0
        DO I=1,NPROCS
           IF( (0 .LT. TMP(I)) .and. (TMP(I) .LE. MGL)) THEN
              NDOM=NDOM+1
              
              IF(NODE_NORTHPOLE_GL.EQ.0) THEN
                 ! Set the north pole global node
                 NODE_NORTHPOLE_GL = TMP(I)
              ELSE
                 ! Make sure it is the same as you global node number
                 IF (NODE_NORTHPOLE_GL.NE.TMP(I)) CALL FATAL_ERROR &
                      &("TWO DOMAINS REPORT DIFFERENT GLOBAL NORTH POLE NODES?")
              END IF
           END IF
        END DO
        DEALLOCATE(TMP)
     END IF
# endif     

     if(dbg_set(dbg_log)) THEN
        WRITE(IPT,*) "! //////////////////////////////////////////////"
        IF(ndom .eq.0) THEN
           WRITE(IPT,*) "! NO NORTH POLE FOUND; PROCEED WITH SPHERICAL MODEL!"
        ELSE
           WRITE(IPT,*) "! FOUND THE GLOBAL NORTH POLE NODE::", NODE_NORTHPOLE_GL
           IF(NDOM .gt. 1)THEN
              WRITE(IPT,*) "THE NORTH POLE IS ON A PROCESSOR BOUNDARY"
              WRITE(IPT,*) "Be afraid, be very afraid...."
              CALL FATAL_ERROR("THE NORTH POLE IS ON A PROCESSOR BOUNDARY")
           END IF
        END IF
        WRITE(IPT,*) "! //////////////////////////////////////////////"
     end if

!     ! ALL OTHER PROCESSORS ESCAPE HERE
!     IF (NODE_NORTHPOLE .EQ. 0) RETURN
     
     ALLOCATE(NODE_NORTHAREA(0:MT)); NODE_NORTHAREA = 0
     ALLOCATE(CELL_NORTHAREA(0:NT)); CELL_NORTHAREA = 0

     ! ALL OTHER PROCESSORS ESCAPE HERE
     IF (NODE_NORTHPOLE .EQ. 0) RETURN
     

     ITMP = NODE_NORTHPOLE
     ALLOCATE(TMP(MT)); TMP = 0     
     MP = 0
     DO I=1,NTSN(ITMP)-1
        MP = MP + 1
        TMP(MP) = NBSN(ITMP,I)
        NODE_NORTHAREA(NBSN(ITMP,I)) = 1
     END DO
     MP = MP + 1
     TMP(MP) = ITMP
     NODE_NORTHAREA(ITMP) = 1
     
     ALLOCATE(MP_LST(MP))
     MP_LST(1:MP) = TMP(1:MP)
     DEALLOCATE(TMP)
     
     ALLOCATE(TMP(NT)); TMP = 0          
     NP = 0 
     DO I=1,NTVE(ITMP)
        NP = NP + 1
        TMP(NP) = NBVE(ITMP,I)
        CELL_NORTHAREA(NBVE(ITMP,I)) = 1
     END DO
     
     ALLOCATE(NP_LST(NP))
     NP_LST(1:NP) = TMP(1:NP)
     DEALLOCATE(TMP)
     
     
     RETURN
     END SUBROUTINE FIND_NORTHPOLE
!==============================================================================|
     
!==============================================================================|

     SUBROUTINE FIND_CELLSIDE
     
     IMPLICIT NONE
     INTEGER  ::  I,IA,IB
     INTEGER, ALLOCATABLE :: TEMP(:)
     
     ! ALL OTHER PROCESSORS ESCAPE HERE
     IF (NODE_NORTHPOLE .EQ. 0) RETURN


    ALLOCATE(TEMP(NE));  TEMP = ZERO
     NPE = 0
     
     DO I=1,NE
       IA = IEC(I,1)
       IB = IEC(I,2)
       IF(CELL_NORTHAREA(IA) == 1 .OR. CELL_NORTHAREA(IB) == 1)THEN
         NPE = NPE + 1
	 TEMP(NPE) = I
       END IF
     END DO
     
     ALLOCATE(NPEDGE_LST(NPE))
     NPEDGE_LST(1:NPE) = TEMP(1:NPE)
     DEALLOCATE(TEMP)
     
     ALLOCATE(TEMP(NCV));  TEMP = ZERO
     NPCV = 0
     
     DO I=1,NCV
       IA = NIEC(I,1)
       IB = NIEC(I,2)
       IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
         NPCV = NPCV + 1
	 TEMP(NPCV) = I
       END IF
     END DO
     
     ALLOCATE(NCEDGE_LST(NPCV))
     NCEDGE_LST(1:NPCV) = TEMP(1:NPCV)
     DEALLOCATE(TEMP)
     
     RETURN
     END SUBROUTINE FIND_CELLSIDE
       	 
!==============================================================================|

!==============================================================================|
   SUBROUTINE ADVAVE_EDGE_XY(XFLUX,YFLUX,IFCETA)

   USE MOD_OBCS
   IMPLICIT NONE
   INTEGER  :: I,J,K,IA,IB,J1,J2,K1,K2,K3,I1,I2,II
   REAL(SP) :: DIJ,ELIJ,XIJ,YIJ,UIJ,VIJ
   REAL(SP) :: COFA1,COFA2,COFA3,COFA4,COFA5,COFA6,COFA7,COFA8
   REAL(SP) :: XADV,YADV,TXXIJ,TYYIJ,TXYIJ 
   REAL(SP) :: VISCOF,VISCOF1,VISCOF2,TEMP
   REAL(SP) :: XFLUX(0:NT),YFLUX(0:NT)
   REAL(SP) :: FACT,FM1,ISWETTMP

   REAL(SP) :: TPA,TPB

   REAL(SP) :: UIJ1_TMP,VIJ1_TMP,UIJ2_TMP,VIJ2_TMP,TXXIJ_TMP,TYYIJ_TMP
   REAL(SP) :: XADV_TMP,YADV_TMP,PSTX_TMP,PSTY_TMP
   REAL(SP) :: UIJ_TMP,VIJ_TMP,UN_TMP
   REAL(SP) :: DLTXC_TMP,DLTYC_TMP
   REAL(SP) :: VX1_TMP,VX2_TMP,VY1_TMP,VY2_TMP
   REAL(SP) :: UAIA,VAIA,UAIB,VAIB,UAK1,VAK1,UAK2,VAK2,UAK3,VAK3
   REAL(SP) :: XIJC_TMP,YIJC_TMP,XCIA_TMP,YCIA_TMP,XCIB_TMP,YCIB_TMP
   REAL(SP) :: XIJ_TMP,YIJ_TMP
   
   REAL(SP) :: IFCETA
  
#  if defined (LIMITED_NO)
   REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
#  else
   REAL(SP),ALLOCATABLE,DIMENSION(:) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
   REAL(SP),ALLOCATABLE,DIMENSION(:) :: UALFA,VALFA
   REAL(SP) :: UALFA_TMP,VALFA_TMP
   INTEGER :: ERROR
   REAL(SP) :: EPS
#  endif  
!------------------------------------------------------------------------------!

   ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN
 
   if(dbg_set(dbg_sbr)) write(ipt,*) "Start: advave_edge_XY"

   SELECT CASE(HORIZONTAL_MIXING_TYPE)
   CASE ('closure')
      FACT = 1.0_SP
      FM1  = 0.0_SP
   CASE('constant')
      FACT = 0.0_SP
      FM1  = 1.0_SP
   CASE DEFAULT
      CALL FATAL_ERROR("UNKNOW HORIZONTAL MIXING TYPE:",&
           & TRIM(HORIZONTAL_MIXING_TYPE) )
   END SELECT

!
!--Initialize Variables--------------------------------------------------------|
!
   DO II=1,NPE
     I=NPEDGE_LST(II)
     IA=IEC(I,1)
     IB=IEC(I,2)
     IF(CELL_NORTHAREA(IA) == 1)THEN
       XFLUX(IA) = 0.0_SP
       YFLUX(IA) = 0.0_SP
       PSTX(IA)  = 0.0_SP
       PSTY(IA)  = 0.0_SP
     END IF  
     IF(CELL_NORTHAREA(IB) == 1)THEN  
       XFLUX(IB) = 0.0_SP
       YFLUX(IB) = 0.0_SP
       PSTX(IB)  = 0.0_SP
       PSTY(IB)  = 0.0_SP
     END IF  
   END DO  
!
!-------------------------ACCUMULATE FLUX OVER ELEMENT EDGES-------------------!
!
#  if !defined (LIMITED_NO)
   ALLOCATE(UIJ1(NPE),VIJ1(NPE),UIJ2(NPE),VIJ2(NPE),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays UIJ1,VIJ1,UIJ2 and VIJ2 can not be allocated", &
   &                  "in SUBROUTINE ADVAVE_XY.")
   
   ALLOCATE(UALFA(0:NT),VALFA(0:NT),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays UALFA,VALFA can not be allocated in SUBROUTINE ADVAVE_XY.")
   
   ALLOCATE(FXX(NPE),FYY(NPE),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays FXX,FYY can not be allocated in SUBROUTINE ADVAVE_XY.")
   
   UIJ1=0.0_SP;VIJ1=0.0_SP;UIJ2=0.0_SP;VIJ2=0.0_SP
   UALFA=1.0_SP;VALFA=1.0_SP
   FXX=0.0_SP;FYY=0.0_SP

   DO II=1,NPE
     I=NPEDGE_LST(II)
     IA=IEC(I,1)
     IB=IEC(I,2)
     J1=IENODE(I,1)
     J2=IENODE(I,2)

#    if !defined (SEMI_IMPLICIT)
     DIJ=0.5_SP*(D(J1)+D(J2))
#    else
     DIJ=0.5_SP*(DT(J1)+DT(J2))
#    endif

!    FLUX FROM LEFT
     K1=NBE(IA,1)
     K2=NBE(IA,2)
     K3=NBE(IA,3)
         
     UAIA = -VA(IA)*COS(XC(IA)*DEG2RAD)-UA(IA)*SIN(XC(IA)*DEG2RAD)
     VAIA = -VA(IA)*SIN(XC(IA)*DEG2RAD)+UA(IA)*COS(XC(IA)*DEG2RAD)
     UAIB = -VA(IB)*COS(XC(IB)*DEG2RAD)-UA(IB)*SIN(XC(IB)*DEG2RAD)
     VAIB = -VA(IB)*SIN(XC(IB)*DEG2RAD)+UA(IB)*COS(XC(IB)*DEG2RAD)

     UAK1 = -VA(K1)*COS(XC(K1)*DEG2RAD)-UA(K1)*SIN(XC(K1)*DEG2RAD)
     VAK1 = -VA(K1)*SIN(XC(K1)*DEG2RAD)+UA(K1)*COS(XC(K1)*DEG2RAD)
     UAK2 = -VA(K2)*COS(XC(K2)*DEG2RAD)-UA(K2)*SIN(XC(K2)*DEG2RAD)
     VAK2 = -VA(K2)*SIN(XC(K2)*DEG2RAD)+UA(K2)*COS(XC(K2)*DEG2RAD)
     UAK3 = -VA(K3)*COS(XC(K3)*DEG2RAD)-UA(K3)*SIN(XC(K3)*DEG2RAD)
     VAK3 = -VA(K3)*SIN(XC(K3)*DEG2RAD)+UA(K3)*COS(XC(K3)*DEG2RAD)
     
     COFA1=A1U_XY(IA,1)*UAIA+A1U_XY(IA,2)*UAK1   &
          +A1U_XY(IA,3)*UAK2+A1U_XY(IA,4)*UAK3
     COFA2=A2U_XY(IA,1)*UAIA+A2U_XY(IA,2)*UAK1   &
          +A2U_XY(IA,3)*UAK2+A2U_XY(IA,4)*UAK3
     COFA5=A1U_XY(IA,1)*VAIA+A1U_XY(IA,2)*VAK1   &
          +A1U_XY(IA,3)*VAK2+A1U_XY(IA,4)*VAK3
     COFA6=A2U_XY(IA,1)*VAIA+A2U_XY(IA,2)*VAK1   &
          +A2U_XY(IA,3)*VAK2+A2U_XY(IA,4)*VAK3
     
     XIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * COS(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
     YIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * SIN(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
     XCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * COS(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
     YCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * SIN(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
 
     XIJ_TMP = XIJC_TMP-XCIA_TMP
     YIJ_TMP = YIJC_TMP-YCIA_TMP

     UIJ1(II)=COFA1*XIJ_TMP+COFA2*YIJ_TMP
     VIJ1(II)=COFA5*XIJ_TMP+COFA6*XIJ_TMP
     UALFA_TMP=ABS(UAIA-UAIB)/ABS(UIJ1(II)+EPSILON(EPS))
     VALFA_TMP=ABS(VAIA-VAIB)/ABS(VIJ1(II)+EPSILON(EPS))
     IF(UALFA_TMP > 1)UALFA_TMP = 1.0_SP
     IF(VALFA_TMP > 1)VALFA_TMP = 1.0_SP
     UALFA(IA)=MIN(UALFA(IA),UALFA_TMP)
     VALFA(IA)=MIN(VALFA(IA),VALFA_TMP)
     
!    FLUX FROM RIGHT
     K1=NBE(IB,1)
     K2=NBE(IB,2)
     K3=NBE(IB,3)
          
!     UAIB = -VA(IB)*COS(XC(IB)*DEG2RAD)-UA(IB)*SIN(XC(IB)*DEG2RAD)
!     VAIB = -VA(IB)*SIN(XC(IB)*DEG2RAD)+UA(IB)*COS(XC(IB)*DEG2RAD)
     UAK1 = -VA(K1)*COS(XC(K1)*DEG2RAD)-UA(K1)*SIN(XC(K1)*DEG2RAD)
     VAK1 = -VA(K1)*SIN(XC(K1)*DEG2RAD)+UA(K1)*COS(XC(K1)*DEG2RAD)
     UAK2 = -VA(K2)*COS(XC(K2)*DEG2RAD)-UA(K2)*SIN(XC(K2)*DEG2RAD)
     VAK2 = -VA(K2)*SIN(XC(K2)*DEG2RAD)+UA(K2)*COS(XC(K2)*DEG2RAD)
     UAK3 = -VA(K3)*COS(XC(K3)*DEG2RAD)-UA(K3)*SIN(XC(K3)*DEG2RAD)
     VAK3 = -VA(K3)*SIN(XC(K3)*DEG2RAD)+UA(K3)*COS(XC(K3)*DEG2RAD)
     
     COFA3=A1U_XY(IB,1)*UAIB+A1U_XY(IB,2)*UAK1   &
          +A1U_XY(IB,3)*UAK2+A1U_XY(IB,4)*UAK3
     COFA4=A2U_XY(IB,1)*UAIB+A2U_XY(IB,2)*UAK1   &
          +A2U_XY(IB,3)*UAK2+A2U_XY(IB,4)*UAK3
     COFA7=A1U_XY(IB,1)*VAIB+A1U_XY(IB,2)*VAK1   &
          +A1U_XY(IB,3)*VAK2+A1U_XY(IB,4)*VAK3
     COFA8=A2U_XY(IB,1)*VAIB+A2U_XY(IB,2)*VAK1   &
          +A2U_XY(IB,3)*VAK2+A2U_XY(IB,4)*VAK3
     
     XCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * COS(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
     YCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * SIN(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
  
     XIJ_TMP = XIJC_TMP-XCIB_TMP
     YIJ_TMP = YIJC_TMP-YCIB_TMP
     UIJ2(II)=COFA3*XIJ_TMP+COFA4*YIJ_TMP
     VIJ2(II)=COFA7*XIJ_TMP+COFA8*YIJ_TMP
     UALFA_TMP=ABS(UAIA-UAIB)/ABS(UIJ2(II)+EPSILON(EPS))
     VALFA_TMP=ABS(VAIA-VAIB)/ABS(VIJ2(II)+EPSILON(EPS))
     IF(UALFA_TMP > 1)UALFA_TMP = 1.0_SP
     IF(VALFA_TMP > 1)VALFA_TMP = 1.0_SP
     UALFA(IB)=MIN(UALFA(IB),UALFA_TMP)
     VALFA(IB)=MIN(VALFA(IB),VALFA_TMP)
     
!    VISCOSITY COEFFICIENT
     VISCOF1=ART(IA)*SQRT(COFA1**2+COFA6**2+0.5_SP*(COFA2+COFA5)**2)
     VISCOF2=ART(IB)*SQRT(COFA3**2+COFA8**2+0.5_SP*(COFA4+COFA7)**2)

!     VISCOF=HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2)/HPRNU + FM1)
!     VISCOF=HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2) + FM1)/HPRNU
     ! David moved HPRNU and added VHC
     VISCOF=(FACT*0.5_SP*(VISCOF1*CC_HVC(IA)+VISCOF2*CC_HVC(IB)) + FM1*0.5_SP*(CC_HVC(IA)+CC_HVC(IB)))/HPRNU


     VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
     VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))

     VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
     VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))

     DLTXC_TMP = VX2_TMP-VX1_TMP
     DLTYC_TMP = VY2_TMP-VY1_TMP

!    SHEAR STRESSES         
     TXXIJ=(COFA1+COFA3)*VISCOF
     TYYIJ=(COFA6+COFA8)*VISCOF
     TXYIJ=0.5_SP*(COFA2+COFA4+COFA5+COFA7)*VISCOF
     FXX(II)=DIJ*(TXXIJ*DLTYC_TMP-TXYIJ*DLTXC_TMP)
     FYY(II)=DIJ*(TXYIJ*DLTYC_TMP-TYYIJ*DLTXC_TMP)
   ENDDO

   DO II=1,NPE
     I=NPEDGE_LST(II)
     IA=IEC(I,1)
     IB=IEC(I,2)
     J1=IENODE(I,1)
     J2=IENODE(I,2)
# if !defined (SEMI_IMPLICIT) 
     DIJ=0.5_SP*(D(J1)+D(J2))
     ELIJ=0.5_SP*(EL(J1)+EL(J2))

!   ggao 0103-2007   consider the sea surface pressure change
# if defined (AIR_PRESSURE)
     ELIJ =ELIJ-0.5_SP*(EL_AIR(J1)+EL_AIR(J2))   !*RAMP
# endif
!     change end

#    if defined (EQUI_TIDE)
     ELIJ=ELIJ-0.5_SP*(EL_EQI(J1)+EL_EQI(J2))
#    endif
#    if defined (ATMO_TIDE)
     ELIJ=ELIJ-0.5_SP*(EL_ATMO(J1)+EL_ATMO(J2))
#    endif

# else

     DIJ= 0.5_SP*(DT(J1)+DT(J2))
     ELIJ=0.5_SP*(ET(J1)+ET(J2))
#    if defined (AIR_PRESSURE)
     ELIJ=ELIJ-( (1.0_SP-IFCETA)*0.5_SP*(EL_AIR(J1)+EL_AIR(J2))+IFCETA*0.5_SP*(ELF_AIR(J1)+ELF_AIR(J2)) )  !*RAMP
#    endif
#    if defined (EQUI_TIDE)
     ELIJ=ELIJ-( (1.0_SP-IFCETA)*0.5_SP*(EL_EQI(J1)+EL_EQI(J2))+IFCETA*0.5_SP*(ELF_EQI(J1)+ELF_EQI(J2)) )
#    endif
#    if defined (ATMO_TIDE)
     ELIJ=ELIJ-( (1.0_SP-IFCETA)*0.5_SP*(EL_ATMO(J1)+EL_ATMO(J2))+IFCETA*0.5_SP*(ELF_ATMO(J1)+ELF_ATMO(J2)) )
#    endif

# endif

     UAIA = -VA(IA)*COS(XC(IA)*DEG2RAD)-UA(IA)*SIN(XC(IA)*DEG2RAD)
     VAIA = -VA(IA)*SIN(XC(IA)*DEG2RAD)+UA(IA)*COS(XC(IA)*DEG2RAD)
     UAIB = -VA(IB)*COS(XC(IB)*DEG2RAD)-UA(IB)*SIN(XC(IB)*DEG2RAD)
     VAIB = -VA(IB)*SIN(XC(IB)*DEG2RAD)+UA(IB)*COS(XC(IB)*DEG2RAD)

     VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
     VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))

     VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
     VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))

     DLTXC_TMP = VX2_TMP-VX1_TMP
     DLTYC_TMP = VY2_TMP-VY1_TMP

     UIJ1(II)=UAIA+UALFA(IA)*UIJ1(II)
     VIJ1(II)=VAIA+VALFA(IA)*VIJ1(II)
     UIJ2(II)=UAIB+UALFA(IB)*UIJ2(II)
     VIJ2(II)=VAIB+VALFA(IB)*VIJ2(II)

#    if defined (LIMITED_1)
     IF(UIJ1(II) > MAX(UAIA,UAIB) .OR. UIJ1(II) < MIN(UAIA,UAIB) .OR.  &
        UIJ2(II) > MAX(UAIA,UAIB) .OR. UIJ2(II) < MIN(UAIA,UAIB))THEN
       UIJ1(II)=UAIA
       UIJ2(II)=UAIB
     END IF

     IF(VIJ1(II) > MAX(VAIA,VAIB) .OR. VIJ1(II) < MIN(VAIA,VAIB) .OR.  &
        VIJ2(II) > MAX(VAIA,VAIB) .OR. VIJ2(II) < MIN(VAIA,VAIB))THEN
       VIJ1(II)=VAIA
       VIJ2(II)=VAIB
     END IF
#    endif
     
     UIJ1_TMP = UIJ1(II)
     VIJ1_TMP = VIJ1(II)
     UIJ2_TMP = UIJ2(II)
     VIJ2_TMP = VIJ2(II)

     UIJ_TMP=0.5_SP*(UIJ1(II)+UIJ2(II))
     VIJ_TMP=0.5_SP*(VIJ1(II)+VIJ2(II))
     UN_TMP=-UIJ_TMP*DLTYC_TMP + VIJ_TMP*DLTXC_TMP

!    ADD CONVECTIVE AND VISCOUS FLUXES

!    ACCUMULATE FLUX
     XADV_TMP=DIJ*UN_TMP*&
              ((1.0_SP-SIGN(1.0_SP,UN_TMP))*UIJ2_TMP                    &
              +(1.0_SP+SIGN(1.0_SP,UN_TMP))*UIJ1_TMP)*0.5_SP  
     YADV_TMP=DIJ*UN_TMP* &
              ((1.0_SP-SIGN(1.0_SP,UN_TMP))*VIJ2_TMP                    &
              +(1.0_SP+SIGN(1.0_SP,UN_TMP))*VIJ1_TMP)*0.5_SP  

     IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
       XFLUX(IA)=XFLUX(IA)+(XADV_TMP+FXX(II)*EPOR(IA))*(1.0_SP-ISBC(I))*IUCP(IA)
       YFLUX(IA)=YFLUX(IA)+(YADV_TMP+FYY(II)*EPOR(IA))*(1.0_SP-ISBC(I))*IUCP(IA)
       XFLUX(IB)=XFLUX(IB)-(XADV_TMP+FXX(II)*EPOR(IB))*(1.0_SP-ISBC(I))*IUCP(IB)
       YFLUX(IB)=YFLUX(IB)-(YADV_TMP+FYY(II)*EPOR(IB))*(1.0_SP-ISBC(I))*IUCP(IB)
     ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
       XFLUX(IA)=XFLUX(IA)+(XADV_TMP+FXX(II)*EPOR(IA))*(1.0_SP-ISBC(I))*IUCP(IA)
       YFLUX(IA)=YFLUX(IA)+(YADV_TMP+FYY(II)*EPOR(IA))*(1.0_SP-ISBC(I))*IUCP(IA)
     ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN
       XFLUX(IB)=XFLUX(IB)-(XADV_TMP+FXX(II)*EPOR(IB))*(1.0_SP-ISBC(I))*IUCP(IB)
       YFLUX(IB)=YFLUX(IB)-(YADV_TMP+FYY(II)*EPOR(IB))*(1.0_SP-ISBC(I))*IUCP(IB)
     END IF 


!    ACCUMULATE BAROTROPIC FLUX

!     VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
!               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
!     VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
!               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))

!     VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
!               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
!     VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
!               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))

!     DLTXC_TMP = VX2_TMP-VX1_TMP
!     DLTYC_TMP = VY2_TMP-VY1_TMP


     IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
#     if !defined (SEMI_IMPLICIT)
      PSTX(IA)=PSTX(IA)-GRAV_E(IA)*D1(IA)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTY(IA)=PSTY(IA)+GRAV_E(IA)*D1(IA)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTX(IB)=PSTX(IB)+GRAV_E(IB)*D1(IB)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
      PSTY(IB)=PSTY(IB)-GRAV_E(IB)*D1(IB)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
#     else
      PSTX(IA)=PSTX(IA)-GRAV_E(IA)*DT1(IA)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTY(IA)=PSTY(IA)+GRAV_E(IA)*DT1(IA)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTX(IB)=PSTX(IB)+GRAV_E(IB)*DT1(IB)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
      PSTY(IB)=PSTY(IB)-GRAV_E(IB)*DT1(IB)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
#     endif
     ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
#     if !defined (SEMI_IMPLICIT)
      PSTX(IA)=PSTX(IA)-GRAV_E(IA)*D1(IA)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTY(IA)=PSTY(IA)+GRAV_E(IA)*D1(IA)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
#     else
      PSTX(IA)=PSTX(IA)-GRAV_E(IA)*DT1(IA)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTY(IA)=PSTY(IA)+GRAV_E(IA)*DT1(IA)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
#     endif
     ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN  
#     if !defined (SEMI_IMPLICIT)
      PSTX(IB)=PSTX(IB)+GRAV_E(IB)*D1(IB)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
      PSTY(IB)=PSTY(IB)-GRAV_E(IB)*D1(IB)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
#     else
      PSTX(IB)=PSTX(IB)+GRAV_E(IB)*DT1(IB)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
      PSTY(IB)=PSTY(IB)-GRAV_E(IB)*DT1(IB)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
#     endif
     END IF

   END DO

   DEALLOCATE(UIJ1,VIJ1,UIJ2,VIJ2,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for UIJ1,VIJ1,UIJ2 and VIJ2", &
   &                  "in subroutine ADVAVE_XY.")
   DEALLOCATE(UALFA,VALFA,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for UALFA,VALFA", &
   &                  "in subroutine ADVAVE_XY.")
   DEALLOCATE(FXX,FYY,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for FXX,FYY", &
   &                  "in subroutine ADVAVE_XY.")
   
#  else

   DO II=1,NPE
     I=NPEDGE_LST(II)
     IA=IEC(I,1)
     IB=IEC(I,2)
     J1=IENODE(I,1)
     J2=IENODE(I,2)

# if !defined (SEMI_IMPLICIT)
     DIJ=0.5_SP*(D(J1)+D(J2))
     ELIJ=0.5_SP*(EL(J1)+EL(J2))

!   ggao 0103-2007   consider the sea surface pressure change
# if defined (AIR_PRESSURE)
     ELIJ =ELIJ-0.5_SP*(EL_AIR(J1)+EL_AIR(J2))   !*RAMP
# endif
!     change end

#    if defined (EQUI_TIDE)
     ELIJ=ELIJ-0.5_SP*(EL_EQI(J1)+EL_EQI(J2))
#    endif
#    if defined (ATMO_TIDE)
     ELIJ=ELIJ-0.5_SP*(EL_ATMO(J1)+EL_ATMO(J2))
#    endif       

# else

     DIJ= 0.5_SP*(DT(J1)+DT(J2))
     ELIJ=0.5_SP*(ET(J1)+ET(J2))
#    if defined (AIR_PRESSURE)
     ELIJ=ELIJ-( (1.0_SP-IFCETA)*0.5_SP*(EL_AIR(J1)+EL_AIR(J2))+IFCETA*0.5_SP*(ELF_AIR(J1)+ELF_AIR(J2)) )  !*RAMP
#    endif
#    if defined (EQUI_TIDE)
     ELIJ=ELIJ-( (1.0_SP-IFCETA)*0.5_SP*(EL_EQI(J1)+EL_EQI(J2))+IFCETA*0.5_SP*(ELF_EQI(J1)+ELF_EQI(J2)) )
#    endif
#    if defined (ATMO_TIDE)
     ELIJ=ELIJ-( (1.0_SP-IFCETA)*0.5_SP*(EL_ATMO(J1)+EL_ATMO(J2))+IFCETA*0.5_SP*(ELF_ATMO(J1)+ELF_ATMO(J2)) )
#    endif     
  
# endif

!    FLUX FROM LEFT
     K1=NBE(IA,1)
     K2=NBE(IA,2)
     K3=NBE(IA,3)
         
     UAIA = -VA(IA)*COS(XC(IA)*DEG2RAD)-UA(IA)*SIN(XC(IA)*DEG2RAD)
     VAIA = -VA(IA)*SIN(XC(IA)*DEG2RAD)+UA(IA)*COS(XC(IA)*DEG2RAD)
     UAK1 = -VA(K1)*COS(XC(K1)*DEG2RAD)-UA(K1)*SIN(XC(K1)*DEG2RAD)
     VAK1 = -VA(K1)*SIN(XC(K1)*DEG2RAD)+UA(K1)*COS(XC(K1)*DEG2RAD)
     UAK2 = -VA(K2)*COS(XC(K2)*DEG2RAD)-UA(K2)*SIN(XC(K2)*DEG2RAD)
     VAK2 = -VA(K2)*SIN(XC(K2)*DEG2RAD)+UA(K2)*COS(XC(K2)*DEG2RAD)
     UAK3 = -VA(K3)*COS(XC(K3)*DEG2RAD)-UA(K3)*SIN(XC(K3)*DEG2RAD)
     VAK3 = -VA(K3)*SIN(XC(K3)*DEG2RAD)+UA(K3)*COS(XC(K3)*DEG2RAD)
     
     COFA1=A1U_XY(IA,1)*UAIA+A1U_XY(IA,2)*UAK1   &
          +A1U_XY(IA,3)*UAK2+A1U_XY(IA,4)*UAK3
     COFA2=A2U_XY(IA,1)*UAIA+A2U_XY(IA,2)*UAK1   &
          +A2U_XY(IA,3)*UAK2+A2U_XY(IA,4)*UAK3
     COFA5=A1U_XY(IA,1)*VAIA+A1U_XY(IA,2)*VAK1   &
          +A1U_XY(IA,3)*VAK2+A1U_XY(IA,4)*VAK3
     COFA6=A2U_XY(IA,1)*VAIA+A2U_XY(IA,2)*VAK1   &
          +A2U_XY(IA,3)*VAK2+A2U_XY(IA,4)*VAK3
     
     XIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * COS(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
     YIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * SIN(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
     XCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * COS(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
     YCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * SIN(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
		  
     XIJ_TMP = XIJC_TMP-XCIA_TMP
     YIJ_TMP = YIJC_TMP-YCIA_TMP

     UIJ1=UAIA+COFA1*XIJ_TMP+COFA2*YIJ_TMP
     VIJ1=VAIA+COFA5*XIJ_TMP+COFA6*XIJ_TMP

!    FLUX FROM RIGHT
     K1=NBE(IB,1)
     K2=NBE(IB,2)
     K3=NBE(IB,3)
          
     UAIB = -VA(IB)*COS(XC(IB)*DEG2RAD)-UA(IB)*SIN(XC(IB)*DEG2RAD)
     VAIB = -VA(IB)*SIN(XC(IB)*DEG2RAD)+UA(IB)*COS(XC(IB)*DEG2RAD)
     UAK1 = -VA(K1)*COS(XC(K1)*DEG2RAD)-UA(K1)*SIN(XC(K1)*DEG2RAD)
     VAK1 = -VA(K1)*SIN(XC(K1)*DEG2RAD)+UA(K1)*COS(XC(K1)*DEG2RAD)
     UAK2 = -VA(K2)*COS(XC(K2)*DEG2RAD)-UA(K2)*SIN(XC(K2)*DEG2RAD)
     VAK2 = -VA(K2)*SIN(XC(K2)*DEG2RAD)+UA(K2)*COS(XC(K2)*DEG2RAD)
     UAK3 = -VA(K3)*COS(XC(K3)*DEG2RAD)-UA(K3)*SIN(XC(K3)*DEG2RAD)
     VAK3 = -VA(K3)*SIN(XC(K3)*DEG2RAD)+UA(K3)*COS(XC(K3)*DEG2RAD)
     
     COFA3=A1U_XY(IB,1)*UAIB+A1U_XY(IB,2)*UAK1   &
          +A1U_XY(IB,3)*UAK2+A1U_XY(IB,4)*UAK3
     COFA4=A2U_XY(IB,1)*UAIB+A2U_XY(IB,2)*UAK1   &
          +A2U_XY(IB,3)*UAK2+A2U_XY(IB,4)*UAK3
     COFA7=A1U_XY(IB,1)*VAIB+A1U_XY(IB,2)*VAK1   &
          +A1U_XY(IB,3)*VAK2+A1U_XY(IB,4)*VAK3
     COFA8=A2U_XY(IB,1)*VAIB+A2U_XY(IB,2)*VAK1   &
          +A2U_XY(IB,3)*VAK2+A2U_XY(IB,4)*VAK3
     
     XCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * COS(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
     YCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * SIN(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
		  
     XIJ_TMP = XIJC_TMP-XCIB_TMP
     YIJ_TMP = YIJC_TMP-YCIB_TMP
     UIJ2=UAIB+COFA3*XIJ_TMP+COFA4*YIJ_TMP
     VIJ2=VAIB+COFA7*XIJ_TMP+COFA8*YIJ_TMP

!    VISCOSITY COEFFICIENT
     VISCOF1=ART(IA)*SQRT(COFA1**2+COFA6**2+0.5_SP*(COFA2+COFA5)**2)
     VISCOF2=ART(IB)*SQRT(COFA3**2+COFA8**2+0.5_SP*(COFA4+COFA7)**2)

     !VISCOF=HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2)/HPRNU + FM1)
     ! David moved HPRNU and added VHC
     VISCOF=(FACT*0.5_SP*(VISCOF1*CC_HVC(IA)+VISCOF2*CC_HVC(IB)) + FM1*0.5_SP*(CC_HVC(IA)+CC_HVC(IB)))/HPRNU
     

     VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
     VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))

     VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
     VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))

     DLTXC_TMP = VX2_TMP-VX1_TMP
     DLTYC_TMP = VY2_TMP-VY1_TMP

!    SHEAR STRESSES         
     TXXIJ=(COFA1+COFA3)*VISCOF
     TYYIJ=(COFA6+COFA8)*VISCOF
     TXYIJ=0.5_SP*(COFA2+COFA4+COFA5+COFA7)*VISCOF
     FXX=DIJ*(TXXIJ*DLTYC_TMP-TXYIJ*DLTXC_TMP)
     FYY=DIJ*(TXYIJ*DLTYC_TMP-TYYIJ*DLTXC_TMP)

     UIJ1_TMP = UIJ1
     VIJ1_TMP = VIJ1
     UIJ2_TMP = UIJ2
     VIJ2_TMP = VIJ2

     UIJ_TMP=0.5_SP*(UIJ1+UIJ2)
     VIJ_TMP=0.5_SP*(VIJ1+VIJ2)
     UN_TMP=-UIJ_TMP*DLTYC_TMP + VIJ_TMP*DLTXC_TMP

!    ADD CONVECTIVE AND VISCOUS FLUXES

!    ACCUMULATE FLUX
     XADV_TMP=DIJ*UN_TMP*&
              ((1.0_SP-SIGN(1.0_SP,UN_TMP))*UIJ2_TMP                    &
              +(1.0_SP+SIGN(1.0_SP,UN_TMP))*UIJ1_TMP)*0.5_SP  
     YADV_TMP=DIJ*UN_TMP* &
              ((1.0_SP-SIGN(1.0_SP,UN_TMP))*VIJ2_TMP                    &
	      +(1.0_SP+SIGN(1.0_SP,UN_TMP))*VIJ1_TMP)*0.5_SP  

     IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
       XFLUX(IA)=XFLUX(IA)+(XADV_TMP+FXX*EPOR(IA))*(1.0_SP-ISBC(I))*IUCP(IA)
       YFLUX(IA)=YFLUX(IA)+(YADV_TMP+FYY*EPOR(IA))*(1.0_SP-ISBC(I))*IUCP(IA)
       XFLUX(IB)=XFLUX(IB)-(XADV_TMP+FXX*EPOR(IB))*(1.0_SP-ISBC(I))*IUCP(IB)
       YFLUX(IB)=YFLUX(IB)-(YADV_TMP+FYY*EPOR(IB))*(1.0_SP-ISBC(I))*IUCP(IB)
     ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
       XFLUX(IA)=XFLUX(IA)+(XADV_TMP+FXX*EPOR(IA))*(1.0_SP-ISBC(I))*IUCP(IA)
       YFLUX(IA)=YFLUX(IA)+(YADV_TMP+FYY*EPOR(IA))*(1.0_SP-ISBC(I))*IUCP(IA)
     ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN
       XFLUX(IB)=XFLUX(IB)-(XADV_TMP+FXX*EPOR(IB))*(1.0_SP-ISBC(I))*IUCP(IB)
       YFLUX(IB)=YFLUX(IB)-(YADV_TMP+FYY*EPOR(IB))*(1.0_SP-ISBC(I))*IUCP(IB)
     END IF 


!    ACCUMULATE BAROTROPIC FLUX

!     VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
!               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
!     VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
!               * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))

!     VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
!               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
!     VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
!               * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))

!     DLTXC_TMP = VX2_TMP-VX1_TMP
!     DLTYC_TMP = VY2_TMP-VY1_TMP


     IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
#     if !defined (SEMI_IMPLICIT)
      PSTX(IA)=PSTX(IA)-GRAV_E(IA)*D1(IA)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTY(IA)=PSTY(IA)+GRAV_E(IA)*D1(IA)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTX(IB)=PSTX(IB)+GRAV_E(IB)*D1(IB)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
      PSTY(IB)=PSTY(IB)-GRAV_E(IB)*D1(IB)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
#     else
      PSTX(IA)=PSTX(IA)-GRAV_E(IA)*DT1(IA)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTY(IA)=PSTY(IA)+GRAV_E(IA)*DT1(IA)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTX(IB)=PSTX(IB)+GRAV_E(IB)*DT1(IB)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
      PSTY(IB)=PSTY(IB)-GRAV_E(IB)*DT1(IB)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
#     endif
     ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
#     if !defined (SEMI_IMPLICIT)
      PSTX(IA)=PSTX(IA)-GRAV_E(IA)*D1(IA)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTY(IA)=PSTY(IA)+GRAV_E(IA)*D1(IA)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
#     else
      PSTX(IA)=PSTX(IA)-GRAV_E(IA)*DT1(IA)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
      PSTY(IA)=PSTY(IA)+GRAV_E(IA)*DT1(IA)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
#     endif
     ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN  
#     if !defined (SEMI_IMPLICIT)
      PSTX(IB)=PSTX(IB)+GRAV_E(IB)*D1(IB)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
      PSTY(IB)=PSTY(IB)-GRAV_E(IB)*D1(IB)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
#     else
      PSTX(IB)=PSTX(IB)+GRAV_E(IB)*DT1(IB)*ELIJ*DLTYC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
      PSTY(IB)=PSTY(IB)-GRAV_E(IB)*DT1(IB)*ELIJ*DLTXC_TMP/(2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
#     endif
     END IF

   END DO
#  endif

   if(dbg_set(dbg_sbr)) write(ipt,*) "End: advave_edge_XY"

   RETURN
   END SUBROUTINE ADVAVE_EDGE_XY
!==============================================================================|

              	 

!==============================================================================|
   SUBROUTINE ADVECTION_EDGE_XY(XFLUX,YFLUX)

   IMPLICIT NONE
!   REAL(SP), INTENT(OUT), DIMENSION(0:NT,KB) :: XFLUX,YFLUX
   REAL(SP) :: XFLUX(0:NT,KB),YFLUX(0:NT,KB)
   REAL(SP) :: DIJ
   REAL(SP) :: COFA1,COFA2,COFA3,COFA4,COFA5,COFA6,COFA7,COFA8
   REAL(SP) :: XADV,YADV,TXXIJ,TYYIJ,TXYIJ,UN
   REAL(SP) :: VISCOF,VISCOF1,VISCOF2,TEMP,TPA,TPB
   REAL(SP) :: XIJA,YIJA,XIJB,YIJB,UIJ,VIJ
   REAL(SP) :: FACT,FM1
   INTEGER  :: I,IA,IB,J1,J2,K1,K2,K3,K4,K5,K6,K,II,J,I1,I2

   REAL(SP) :: UIJ1_TMP,VIJ1_TMP,UIJ2_TMP,VIJ2_TMP,TXXIJ_TMP,TYYIJ_TMP
   REAL(SP) :: XADV_TMP,YADV_TMP
   REAL(SP) :: UIJ_TMP,VIJ_TMP,UN_TMP
   REAL(SP) :: DLTXC_TMP,DLTYC_TMP
   REAL(SP) :: VX1_TMP,VX2_TMP,VY1_TMP,VY2_TMP
   REAL(SP) :: UIA,VIA,UIB,VIB,UK1,VK1,UK2,VK2,UK3,VK3,UK4,VK4,UK5,VK5,UK6,VK6
   REAL(SP) :: XIJC_TMP,YIJC_TMP,XCIA_TMP,YCIA_TMP,XCIB_TMP,YCIB_TMP

#  if defined (LIMITED_NO)
   REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
#  else
   REAL(SP),ALLOCATABLE,DIMENSION(:) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
   REAL(SP),ALLOCATABLE,DIMENSION(:) :: UALFA,VALFA
   REAL(SP) :: UALFA_TMP,VALFA_TMP
   INTEGER :: ERROR
   REAL(SP) :: EPS
#  endif
!------------------------------------------------------------------------------|
   ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN

   if(dbg_set(dbg_sbr)) write(ipt,*) "Start: advection_edge_XY"

   SELECT CASE(HORIZONTAL_MIXING_TYPE)
   CASE ('closure')
      FACT = 1.0_SP
      FM1  = 0.0_SP
   CASE('constant')
      FACT = 0.0_SP
      FM1  = 1.0_SP
   CASE DEFAULT
      CALL FATAL_ERROR("UNKNOW HORIZONTAL MIXING TYPE:",&
           & TRIM(HORIZONTAL_MIXING_TYPE) )
   END SELECT

!
!--Initialize Variables--------------------------------------------------------|
!
   DO K = 1,KBM1
     DO II=1,NPE
       I=NPEDGE_LST(II)
       IA=IEC(I,1)
       IB=IEC(I,2)
       IF(CELL_NORTHAREA(IA) == 1)THEN
         XFLUX(IA,K) = 0.0_SP
         YFLUX(IA,K) = 0.0_SP
       END IF
       IF(CELL_NORTHAREA(IB) == 1)THEN	 
         XFLUX(IB,K) = 0.0_SP
         YFLUX(IB,K) = 0.0_SP
       END IF	 
     END DO
   END DO    

!
!--Loop Over Edges and Accumulate Fluxes-For Each Element----------------------|
!

#  if !defined (LIMITED_NO)
   ALLOCATE(UIJ1(NPE),VIJ1(NPE),UIJ2(NPE),VIJ2(NPE),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays UIJ1,VIJ1,UIJ2 and VIJ2 can not be allocated", &
   &                  "in SUBROUTINE ADVECTION_XY.")
   
   ALLOCATE(UALFA(0:NT),VALFA(0:NT),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays UALFA,VALFA can not be allocated", &
   &                  "in SUBROUTINE ADVECTION_XY.")
   
   ALLOCATE(FXX(NPE),FYY(NPE),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays FXX,FYY can not be allocated", &
   &                  "in SUBROUTINE ADVECTION_XY.")
   
   DO K=1,KBM1
     UIJ1=0.0_SP;VIJ1=0.0_SP;UIJ2=0.0_SP;VIJ2=0.0_SP
     UALFA=1.0_SP;VALFA=1.0_SP
     FXX=0.0_SP;FYY=0.0_SP
     
     DO II=1,NPE
       I=NPEDGE_LST(II)
       IA=IEC(I,1)
       IB=IEC(I,2)
       J1=IENODE(I,1)
       J2=IENODE(I,2)
!     DIJ= 0.5_SP*(DT(J1)+DT(J2))

       K1=NBE(IA,1)
       K2=NBE(IA,2)
       K3=NBE(IA,3)
       K4=NBE(IB,1)
       K5=NBE(IB,2)
       K6=NBE(IB,3)

       XIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * COS(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
       YIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * SIN(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
       XCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * COS(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
       YCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * SIN(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
       XCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * COS(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
       YCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * SIN(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
  
       XIJA = XIJC_TMP-XCIA_TMP
       YIJA = YIJC_TMP-YCIA_TMP
       XIJB = XIJC_TMP-XCIB_TMP
       YIJB = YIJC_TMP-YCIB_TMP

       DIJ= 0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))

       UIA = -V(IA,K)*COS(XC(IA)*DEG2RAD)-U(IA,K)*SIN(XC(IA)*DEG2RAD)
       VIA = -V(IA,K)*SIN(XC(IA)*DEG2RAD)+U(IA,K)*COS(XC(IA)*DEG2RAD)
       UIB = -V(IB,K)*COS(XC(IB)*DEG2RAD)-U(IB,K)*SIN(XC(IB)*DEG2RAD)
       VIB = -V(IB,K)*SIN(XC(IB)*DEG2RAD)+U(IB,K)*COS(XC(IB)*DEG2RAD)
       UK1 = -V(K1,K)*COS(XC(K1)*DEG2RAD)-U(K1,K)*SIN(XC(K1)*DEG2RAD)
       VK1 = -V(K1,K)*SIN(XC(K1)*DEG2RAD)+U(K1,K)*COS(XC(K1)*DEG2RAD)
       UK2 = -V(K2,K)*COS(XC(K2)*DEG2RAD)-U(K2,K)*SIN(XC(K2)*DEG2RAD)
       VK2 = -V(K2,K)*SIN(XC(K2)*DEG2RAD)+U(K2,K)*COS(XC(K2)*DEG2RAD)
       UK3 = -V(K3,K)*COS(XC(K3)*DEG2RAD)-U(K3,K)*SIN(XC(K3)*DEG2RAD)
       VK3 = -V(K3,K)*SIN(XC(K3)*DEG2RAD)+U(K3,K)*COS(XC(K3)*DEG2RAD)
       UK4 = -V(K4,K)*COS(XC(K4)*DEG2RAD)-U(K4,K)*SIN(XC(K4)*DEG2RAD)
       VK4 = -V(K4,K)*SIN(XC(K4)*DEG2RAD)+U(K4,K)*COS(XC(K4)*DEG2RAD)
       UK5 = -V(K5,K)*COS(XC(K5)*DEG2RAD)-U(K5,K)*SIN(XC(K5)*DEG2RAD)
       VK5 = -V(K5,K)*SIN(XC(K5)*DEG2RAD)+U(K5,K)*COS(XC(K5)*DEG2RAD)
       UK6 = -V(K6,K)*COS(XC(K6)*DEG2RAD)-U(K6,K)*SIN(XC(K6)*DEG2RAD)
       VK6 = -V(K6,K)*SIN(XC(K6)*DEG2RAD)+U(K6,K)*COS(XC(K6)*DEG2RAD)
       !!FORM THE LEFT FLUX
       COFA1=A1U_XY(IA,1)*UIA+A1U_XY(IA,2)*UK1   &
            +A1U_XY(IA,3)*UK2+A1U_XY(IA,4)*UK3
       COFA2=A2U_XY(IA,1)*UIA+A2U_XY(IA,2)*UK1   &
            +A2U_XY(IA,3)*UK2+A2U_XY(IA,4)*UK3
       COFA5=A1U_XY(IA,1)*VIA+A1U_XY(IA,2)*VK1   &
            +A1U_XY(IA,3)*VK2+A1U_XY(IA,4)*VK3
       COFA6=A2U_XY(IA,1)*VIA+A2U_XY(IA,2)*VK1   &
            +A2U_XY(IA,3)*VK2+A2U_XY(IA,4)*VK3

       UIJ1(II)=COFA1*XIJA+COFA2*YIJA
       VIJ1(II)=COFA5*XIJA+COFA6*YIJA
       UALFA_TMP=ABS(UIA-UIB)/ABS(UIJ1(II)+EPSILON(EPS))
       VALFA_TMP=ABS(VIA-VIB)/ABS(VIJ1(II)+EPSILON(EPS))
       IF(UALFA_TMP > 1)UALFA_TMP = 1.0_SP
       IF(VALFA_TMP > 1)VALFA_TMP = 1.0_SP
       UALFA(IA)=MIN(UALFA(IA),UALFA_TMP)
       VALFA(IA)=MIN(VALFA(IA),VALFA_TMP)
       
       !!FORM THE RIGHT FLUX
       COFA3=A1U_XY(IB,1)*UIB+A1U_XY(IB,2)*UK4   &
            +A1U_XY(IB,3)*UK5+A1U_XY(IB,4)*UK6
       COFA4=A2U_XY(IB,1)*UIB+A2U_XY(IB,2)*UK4   &
            +A2U_XY(IB,3)*UK5+A2U_XY(IB,4)*UK6
       COFA7=A1U_XY(IB,1)*VIB+A1U_XY(IB,2)*VK4   &
            +A1U_XY(IB,3)*VK5+A1U_XY(IB,4)*VK6
       COFA8=A2U_XY(IB,1)*VIB+A2U_XY(IB,2)*VK4   &
            +A2U_XY(IB,3)*VK5+A2U_XY(IB,4)*VK6

       UIJ2(II)=COFA3*XIJB+COFA4*YIJB
       VIJ2(II)=COFA7*XIJB+COFA8*YIJB
       UALFA_TMP=ABS(UIA-UIB)/ABS(UIJ2(II)+EPSILON(EPS))
       VALFA_TMP=ABS(VIA-VIB)/ABS(VIJ2(II)+EPSILON(EPS))
       IF(UALFA_TMP > 1)UALFA_TMP = 1.0_SP
       IF(VALFA_TMP > 1)VALFA_TMP = 1.0_SP
       UALFA(IB)=MIN(UALFA(IB),UALFA_TMP)
       VALFA(IB)=MIN(VALFA(IB),VALFA_TMP)

       !!    VISCOSITY COEFFICIENT EDGE
       VISCOF1=ART(IA)*SQRT(COFA1**2+COFA6**2+0.5_SP*(COFA2+COFA5)**2)
       VISCOF2=ART(IB)*SQRT(COFA3**2+COFA8**2+0.5_SP*(COFA4+COFA7)**2)
        
!       VISCOF = HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2)/HPRNU + FM1)
!       VISCOF = HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2) + FM1)/HPRNU
       ! David moved HPRNU and added VHC
       VISCOF=(FACT*0.5_SP*(VISCOF1*CC_HVC(IA)+VISCOF2*CC_HVC(IB)) + FM1*0.5_SP*(CC_HVC(IA)+CC_HVC(IB)))/HPRNU

       VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD)       &
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))
       VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))

       VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))
       VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))

       DLTXC_TMP = VX2_TMP-VX1_TMP
       DLTYC_TMP = VY2_TMP-VY1_TMP
       
       TXXIJ=(COFA1+COFA3)*VISCOF
       TYYIJ=(COFA6+COFA8)*VISCOF
       TXYIJ=0.5_SP*(COFA2+COFA4+COFA5+COFA7)*VISCOF
       FXX(II)=DIJ*(TXXIJ*DLTYC_TMP-TXYIJ*DLTXC_TMP)
       FYY(II)=DIJ*(TXYIJ*DLTYC_TMP-TYYIJ*DLTXC_TMP)
     ENDDO

     DO II=1,NPE
       I=NPEDGE_LST(II)
       IA=IEC(I,1)
       IB=IEC(I,2)
       J1=IENODE(I,1)
       J2=IENODE(I,2)

       DIJ= 0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))
       UIA = -V(IA,K)*COS(XC(IA)*DEG2RAD)-U(IA,K)*SIN(XC(IA)*DEG2RAD)
       VIA = -V(IA,K)*SIN(XC(IA)*DEG2RAD)+U(IA,K)*COS(XC(IA)*DEG2RAD)
       UIB = -V(IB,K)*COS(XC(IB)*DEG2RAD)-U(IB,K)*SIN(XC(IB)*DEG2RAD)
       VIB = -V(IB,K)*SIN(XC(IB)*DEG2RAD)+U(IB,K)*COS(XC(IB)*DEG2RAD)

       VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD)       &
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))
       VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))

       VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))
       VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))

       DLTXC_TMP = VX2_TMP-VX1_TMP
       DLTYC_TMP = VY2_TMP-VY1_TMP

       UIJ1(II)=UIA+UALFA(IA)*UIJ1(II)
       VIJ1(II)=VIA+VALFA(IA)*VIJ1(II)
       UIJ2(II)=UIB+UALFA(IB)*UIJ2(II)
       VIJ2(II)=VIB+VALFA(IB)*VIJ2(II)

#      if defined (LIMITED_1)
       IF(UIJ1(II) > MAX(UIA,UIB) .OR. UIJ1(II) < MIN(UIA,UIB) .OR. &
          UIJ2(II) > MAX(UIA,UIB) .OR. UIJ2(II) < MIN(UIA,UIB))THEN
         UIJ1(II)=UIA
         UIJ2(II)=UIB
       END IF

       IF(VIJ1(II) > MAX(VIA,VIB) .OR. VIJ1(II) < MIN(VIA,VIB) .OR. &
          VIJ2(II) > MAX(VIA,VIB) .OR. VIJ2(II) < MIN(VIA,VIB))THEN
         VIJ1(II)=VIA
         VIJ2(II)=VIB
       END IF
#      endif       

       UIJ1_TMP = UIJ1(II)
       VIJ1_TMP = VIJ1(II)
       UIJ2_TMP = UIJ2(II)
       VIJ2_TMP = VIJ2(II)
       
       UIJ_TMP=0.5_SP*(UIJ1(II)+UIJ2(II))
       VIJ_TMP=0.5_SP*(VIJ1(II)+VIJ2(II))
       UN_TMP=-UIJ_TMP*DLTYC_TMP + VIJ_TMP*DLTXC_TMP

       !!COMPUTE BOUNDARY FLUX AUGMENTERS   
       TPA = FLOAT(1-ISBC(I))*EPOR(IA)
       TPB = FLOAT(1-ISBC(I))*EPOR(IB)


       !!ACCUMULATE THE FLUX
       XADV_TMP=DIJ*UN_TMP*&
                ((1.0_SP-SIGN(1.0_SP,UN_TMP))*UIJ2_TMP                 &
                +(1.0_SP+SIGN(1.0_SP,UN_TMP))*UIJ1_TMP)*0.5_SP
       YADV_TMP=DIJ*UN_TMP* &
                ((1.0_SP-SIGN(1.0_SP,UN_TMP))*VIJ2_TMP                 &
                +(1.0_SP+SIGN(1.0_SP,UN_TMP))*VIJ1_TMP)*0.5_SP 
       IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
         XFLUX(IA,K)=XFLUX(IA,K)+XADV_TMP*TPA+(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX(IA,K)=YFLUX(IA,K)+YADV_TMP*TPA+(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IA)
         XFLUX(IB,K)=XFLUX(IB,K)-XADV_TMP*TPB-(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX(IB,K)=YFLUX(IB,K)-YADV_TMP*TPB-(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IB)
       ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
         XFLUX(IA,K)=XFLUX(IA,K)+XADV_TMP*TPA+(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX(IA,K)=YFLUX(IA,K)+YADV_TMP*TPA+(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IA)
       ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN 
         XFLUX(IB,K)=XFLUX(IB,K)-XADV_TMP*TPB-(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX(IB,K)=YFLUX(IB,K)-YADV_TMP*TPB-(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IB)
       END IF
     END DO
   END DO

   DEALLOCATE(UIJ1,VIJ1,UIJ2,VIJ2,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for UIJ1,VIJ1,UIJ2 and VIJ2", &
   &                  "in subroutine ADVECTION_XY.")
   DEALLOCATE(UALFA,VALFA,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for UALFA,VALFA", &
   &                  "in subroutine ADVECTION_XY.")
   DEALLOCATE(FXX,FYY,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for FXX,FYY", &
   &                  "in subroutine ADVECTION_XY.")

#  else

   DO II=1,NPE
     I=NPEDGE_LST(II)
     IA=IEC(I,1)
     IB=IEC(I,2)
     J1=IENODE(I,1)
     J2=IENODE(I,2)
!     DIJ= 0.5_SP*(DT(J1)+DT(J2))

     K1=NBE(IA,1)
     K2=NBE(IA,2)
     K3=NBE(IA,3)
     K4=NBE(IB,1)
     K5=NBE(IB,2)
     K6=NBE(IB,3)

     XIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * COS(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
     YIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * SIN(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
     XCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * COS(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
     YCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * SIN(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
     XCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * COS(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
     YCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * SIN(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
		  
     XIJA = XIJC_TMP-XCIA_TMP
     YIJA = YIJC_TMP-YCIA_TMP
     XIJB = XIJC_TMP-XCIB_TMP
     YIJB = YIJC_TMP-YCIB_TMP

     DO K=1,KBM1

       DIJ= 0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))

       UIA = -V(IA,K)*COS(XC(IA)*DEG2RAD)-U(IA,K)*SIN(XC(IA)*DEG2RAD)
       VIA = -V(IA,K)*SIN(XC(IA)*DEG2RAD)+U(IA,K)*COS(XC(IA)*DEG2RAD)
       UIB = -V(IB,K)*COS(XC(IB)*DEG2RAD)-U(IB,K)*SIN(XC(IB)*DEG2RAD)
       VIB = -V(IB,K)*SIN(XC(IB)*DEG2RAD)+U(IB,K)*COS(XC(IB)*DEG2RAD)
       UK1 = -V(K1,K)*COS(XC(K1)*DEG2RAD)-U(K1,K)*SIN(XC(K1)*DEG2RAD)
       VK1 = -V(K1,K)*SIN(XC(K1)*DEG2RAD)+U(K1,K)*COS(XC(K1)*DEG2RAD)
       UK2 = -V(K2,K)*COS(XC(K2)*DEG2RAD)-U(K2,K)*SIN(XC(K2)*DEG2RAD)
       VK2 = -V(K2,K)*SIN(XC(K2)*DEG2RAD)+U(K2,K)*COS(XC(K2)*DEG2RAD)
       UK3 = -V(K3,K)*COS(XC(K3)*DEG2RAD)-U(K3,K)*SIN(XC(K3)*DEG2RAD)
       VK3 = -V(K3,K)*SIN(XC(K3)*DEG2RAD)+U(K3,K)*COS(XC(K3)*DEG2RAD)
       UK4 = -V(K4,K)*COS(XC(K4)*DEG2RAD)-U(K4,K)*SIN(XC(K4)*DEG2RAD)
       VK4 = -V(K4,K)*SIN(XC(K4)*DEG2RAD)+U(K4,K)*COS(XC(K4)*DEG2RAD)
       UK5 = -V(K5,K)*COS(XC(K5)*DEG2RAD)-U(K5,K)*SIN(XC(K5)*DEG2RAD)
       VK5 = -V(K5,K)*SIN(XC(K5)*DEG2RAD)+U(K5,K)*COS(XC(K5)*DEG2RAD)
       UK6 = -V(K6,K)*COS(XC(K6)*DEG2RAD)-U(K6,K)*SIN(XC(K6)*DEG2RAD)
       VK6 = -V(K6,K)*SIN(XC(K6)*DEG2RAD)+U(K6,K)*COS(XC(K6)*DEG2RAD)
       !!FORM THE LEFT FLUX
       COFA1=A1U_XY(IA,1)*UIA+A1U_XY(IA,2)*UK1   &
            +A1U_XY(IA,3)*UK2+A1U_XY(IA,4)*UK3
       COFA2=A2U_XY(IA,1)*UIA+A2U_XY(IA,2)*UK1   &
            +A2U_XY(IA,3)*UK2+A2U_XY(IA,4)*UK3
       COFA5=A1U_XY(IA,1)*VIA+A1U_XY(IA,2)*VK1   &
            +A1U_XY(IA,3)*VK2+A1U_XY(IA,4)*VK3
       COFA6=A2U_XY(IA,1)*VIA+A2U_XY(IA,2)*VK1   &
            +A2U_XY(IA,3)*VK2+A2U_XY(IA,4)*VK3
       UIJ1=UIA+COFA1*XIJA+COFA2*YIJA
       VIJ1=VIA+COFA5*XIJA+COFA6*YIJA

       !!FORM THE RIGHT FLUX
       COFA3=A1U_XY(IB,1)*UIB+A1U_XY(IB,2)*UK4   &
            +A1U_XY(IB,3)*UK5+A1U_XY(IB,4)*UK6
       COFA4=A2U_XY(IB,1)*UIB+A2U_XY(IB,2)*UK4   &
            +A2U_XY(IB,3)*UK5+A2U_XY(IB,4)*UK6
       COFA7=A1U_XY(IB,1)*VIB+A1U_XY(IB,2)*VK4   &
            +A1U_XY(IB,3)*VK5+A1U_XY(IB,4)*VK6
       COFA8=A2U_XY(IB,1)*VIB+A2U_XY(IB,2)*VK4   &
            +A2U_XY(IB,3)*VK5+A2U_XY(IB,4)*VK6
       UIJ2=UIB+COFA3*XIJB+COFA4*YIJB
       VIJ2=VIB+COFA7*XIJB+COFA8*YIJB

       !!    VISCOSITY COEFFICIENT EDGE
       VISCOF1=ART(IA)*SQRT(COFA1**2+COFA6**2+0.5_SP*(COFA2+COFA5)**2)
       VISCOF2=ART(IB)*SQRT(COFA3**2+COFA8**2+0.5_SP*(COFA4+COFA7)**2)
        
!       VISCOF = HORCON*(FACT*0.5_SP*(VISCOF1+VISCOF2)/HPRNU + FM1)
       ! David moved HPRNU and added VHC
       VISCOF=(FACT*0.5_SP*(VISCOF1*CC_HVC(IA)+VISCOF2*CC_HVC(IB)) + FM1*0.5_SP*(CC_HVC(IA)+CC_HVC(IB)))/HPRNU

       VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD)       &
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))
       VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))

       VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))
       VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))

       DLTXC_TMP = VX2_TMP-VX1_TMP
       DLTYC_TMP = VY2_TMP-VY1_TMP
       
       TXXIJ=(COFA1+COFA3)*VISCOF
       TYYIJ=(COFA6+COFA8)*VISCOF
       TXYIJ=0.5_SP*(COFA2+COFA4+COFA5+COFA7)*VISCOF
       FXX=DIJ*(TXXIJ*DLTYC_TMP-TXYIJ*DLTXC_TMP)
       FYY=DIJ*(TXYIJ*DLTYC_TMP-TYYIJ*DLTXC_TMP)

       UIJ1_TMP = UIJ1
       VIJ1_TMP = VIJ1
       UIJ2_TMP = UIJ2
       VIJ2_TMP = VIJ2
     
       UIJ_TMP=0.5_SP*(UIJ1+UIJ2)
       VIJ_TMP=0.5_SP*(VIJ1+VIJ2)
       UN_TMP=-UIJ_TMP*DLTYC_TMP + VIJ_TMP*DLTXC_TMP

       !!COMPUTE BOUNDARY FLUX AUGMENTERS   
       TPA = FLOAT(1-ISBC(I))*EPOR(IA)
       TPB = FLOAT(1-ISBC(I))*EPOR(IB)


       !!ACCUMULATE THE FLUX
       XADV_TMP=DIJ*UN_TMP*&
                ((1.0_SP-SIGN(1.0_SP,UN_TMP))*UIJ2_TMP                 &
	        +(1.0_SP+SIGN(1.0_SP,UN_TMP))*UIJ1_TMP)*0.5_SP
       YADV_TMP=DIJ*UN_TMP* &
                ((1.0_SP-SIGN(1.0_SP,UN_TMP))*VIJ2_TMP                 &
	          +(1.0_SP+SIGN(1.0_SP,UN_TMP))*VIJ1_TMP)*0.5_SP 
       IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
         XFLUX(IA,K)=XFLUX(IA,K)+XADV_TMP*TPA+(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX(IA,K)=YFLUX(IA,K)+YADV_TMP*TPA+(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IA)
         XFLUX(IB,K)=XFLUX(IB,K)-XADV_TMP*TPB-(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX(IB,K)=YFLUX(IB,K)-YADV_TMP*TPB-(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IB)
       ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
         XFLUX(IA,K)=XFLUX(IA,K)+XADV_TMP*TPA+(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX(IA,K)=YFLUX(IA,K)+YADV_TMP*TPA+(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IA)
       ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN 
         XFLUX(IB,K)=XFLUX(IB,K)-XADV_TMP*TPB-(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX(IB,K)=YFLUX(IB,K)-YADV_TMP*TPB-(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IB)
       END IF
     END DO
   END DO
#  endif

   if(dbg_set(dbg_sbr)) write(ipt,*) "End: advection_edge_XY"

   RETURN
   END SUBROUTINE ADVECTION_EDGE_XY
!==============================================================================|

!==============================================================================!

   SUBROUTINE ADV_UV_EDGE_XY(XFLUX,YFLUX,CETA,STG,K_STG)

   IMPLICIT NONE
   REAL(SP) :: XFLUX(0:NT,KB),YFLUX(0:NT,KB)
   REAL(SP) :: PSTX_TM(0:NT,KB),PSTY_TM(0:NT,KB)
   REAL(SP) :: COFA1,COFA2,COFA3,COFA4,COFA5,COFA6,COFA7,COFA8
   REAL(SP) :: XADV,YADV,TXXIJ,TYYIJ,TXYIJ,UN
   REAL(SP) :: VISCOF,VISCOF1,VISCOF2,TEMP,TPA,TPB
   REAL(SP) :: XIJA,YIJA,XIJB,YIJB,UIJ,VIJ
   REAL(SP) :: SITA,DIJ,ELIJ,TMPA,TMPB,TMP,XFLUXV,YFLUXV
   REAL(SP) :: FACT,FM1,EXFLUX,ISWETTMP
   INTEGER  :: I,IA,IB,J1,J2,K1,K2,K3,K4,K5,K6,K,II,J,I1,I2

!   REAL(SP) :: TY
   REAL(SP) :: UIJ1_TMP,VIJ1_TMP,UIJ2_TMP,VIJ2_TMP,UIJ3_TMP,VIJ3_TMP
   REAL(SP) :: U_TMP,V_TMP,UF_TMP,VF_TMP
   REAL(SP) :: TXXIJ_TMP,TYYIJ_TMP
   REAL(SP) :: XADV_TMP,YADV_TMP,PSTX_TMP,PSTY_TMP
   REAL(SP) :: UIJ_TMP,VIJ_TMP,EXFLUX_TMP
   REAL(SP) :: DLTXC_TMP,DLTYC_TMP
   REAL(SP) :: VX1_TMP,VX2_TMP,VY1_TMP,VY2_TMP
   REAL(SP) :: UIA,VIA,UIB,VIB,UK1,VK1,UK2,VK2,UK3,VK3,UK4,VK4,UK5,VK5,UK6,VK6 
   REAL(SP) :: XIJC_TMP,YIJC_TMP,XCIA_TMP,YCIA_TMP,XCIB_TMP,YCIB_TMP

   REAL(SP) :: CETA

#  if defined (LIMITED_NO)
   REAL(SP) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
#  else
   REAL(SP),ALLOCATABLE,DIMENSION(:) :: UIJ1,VIJ1,UIJ2,VIJ2,FXX,FYY
   REAL(SP),ALLOCATABLE,DIMENSION(:) :: UALFA,VALFA
   REAL(SP) :: UALFA_TMP,VALFA_TMP
   INTEGER :: ERROR
   REAL(SP) :: EPS
#  endif

   INTEGER :: STG, K_STG
!------------------------------------------------------------------------------!
  ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN
 
   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: adv_uv_edge_xy"

   SELECT CASE(HORIZONTAL_MIXING_TYPE)
   CASE ('closure')
      FACT = 1.0_SP
      FM1  = 0.0_SP
   CASE('constant')
      FACT = 0.0_SP
      FM1  = 1.0_SP
   CASE DEFAULT
      CALL FATAL_ERROR("UNKNOW HORIZONTAL MIXING TYPE:",&
           & TRIM(HORIZONTAL_MIXING_TYPE) )
   END SELECT

!
!-----Initialize Flux Variables------------------------------------------------!
!

   DO K = 1,KBM1
     DO II=1,NPE
       I=NPEDGE_LST(II)
       IA=IEC(I,1)
       IB=IEC(I,2)
       IF(CELL_NORTHAREA(IA) == 1)THEN
         XFLUX(IA,K) = 0.0_SP
         YFLUX(IA,K) = 0.0_SP
         PSTX_TM(IA,K) = 0.0_SP
         PSTY_TM(IA,K) = 0.0_SP
#        if defined (SEMI_IMPLICIT)
         XFLUX3_NP(IA,K) = 0.0_SP
         YFLUX3_NP(IA,K) = 0.0_SP
#        endif
       END IF
       IF(CELL_NORTHAREA(IB) == 1)THEN
         XFLUX(IB,K) = 0.0_SP
         YFLUX(IB,K) = 0.0_SP
         PSTX_TM(IB,K) = 0.0_SP
         PSTY_TM(IB,K) = 0.0_SP
#        if defined (SEMI_IMPLICIT)
         XFLUX3_NP(IB,K) = 0.0_SP
         YFLUX3_NP(IB,K) = 0.0_SP
#        endif
       END IF 
     END DO
   END DO         

!
!-----Loop Over Edges and Accumulate Flux--------------------------------------!
!
#  if !defined (LIMITED_NO)
   ALLOCATE(UIJ1(NPE),VIJ1(NPE),UIJ2(NPE),VIJ2(NPE),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays UIJ1,VIJ1,UIJ2 and VIJ2 can not be allocated", &
   &                  "in SUBROUTINE ADV_UV_EDGE_XY.")
   
   ALLOCATE(UALFA(0:NT),VALFA(0:NT),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays UALFA,VALFA can not be allocated", &
   &                  "in SUBROUTINE ADV_UV_EDGE_XY.")
   
   ALLOCATE(FXX(NPE),FYY(NPE),STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("The arrays FXX,FYY can not be allocated", &
   &                  "in SUBROUTINE ADV_UV_EDGE_XY.")
   
   DO K=1,KBM1
     UIJ1=0.0_SP;VIJ1=0.0_SP;UIJ2=0.0_SP;VIJ2=0.0_SP
     UALFA=1.0_SP;VALFA=1.0_SP
     FXX=0.0_SP;FYY=0.0_SP

     DO II=1,NPE
       I=NPEDGE_LST(II)
       IA=IEC(I,1)
       IB=IEC(I,2)
       J1=IENODE(I,1)
       J2=IENODE(I,2)
!       DIJ=0.5_SP*(DT(J1)+DT(J2))

       K1=NBE(IA,1)
       K2=NBE(IA,2)
       K3=NBE(IA,3)
       K4=NBE(IB,1)
       K5=NBE(IB,2)
       K6=NBE(IB,3)

       XIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * COS(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
       YIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * SIN(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
  
       XCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * COS(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
       YCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * SIN(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
       XCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * COS(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
       YCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * SIN(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
  
       XIJA = XIJC_TMP-XCIA_TMP
       YIJA = YIJC_TMP-YCIA_TMP
       XIJB = XIJC_TMP-XCIB_TMP
       YIJB = YIJC_TMP-YCIB_TMP

       DIJ=0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))

       UIA = -V(IA,K)*COS(XC(IA)*DEG2RAD)-U(IA,K)*SIN(XC(IA)*DEG2RAD)
       VIA = -V(IA,K)*SIN(XC(IA)*DEG2RAD)+U(IA,K)*COS(XC(IA)*DEG2RAD)
       UIB = -V(IB,K)*COS(XC(IB)*DEG2RAD)-U(IB,K)*SIN(XC(IB)*DEG2RAD)
       VIB = -V(IB,K)*SIN(XC(IB)*DEG2RAD)+U(IB,K)*COS(XC(IB)*DEG2RAD)
       UK1 = -V(K1,K)*COS(XC(K1)*DEG2RAD)-U(K1,K)*SIN(XC(K1)*DEG2RAD)
       VK1 = -V(K1,K)*SIN(XC(K1)*DEG2RAD)+U(K1,K)*COS(XC(K1)*DEG2RAD)
       UK2 = -V(K2,K)*COS(XC(K2)*DEG2RAD)-U(K2,K)*SIN(XC(K2)*DEG2RAD)
       VK2 = -V(K2,K)*SIN(XC(K2)*DEG2RAD)+U(K2,K)*COS(XC(K2)*DEG2RAD)
       UK3 = -V(K3,K)*COS(XC(K3)*DEG2RAD)-U(K3,K)*SIN(XC(K3)*DEG2RAD)
       VK3 = -V(K3,K)*SIN(XC(K3)*DEG2RAD)+U(K3,K)*COS(XC(K3)*DEG2RAD)
       UK4 = -V(K4,K)*COS(XC(K4)*DEG2RAD)-U(K4,K)*SIN(XC(K4)*DEG2RAD)
       VK4 = -V(K4,K)*SIN(XC(K4)*DEG2RAD)+U(K4,K)*COS(XC(K4)*DEG2RAD)
       UK5 = -V(K5,K)*COS(XC(K5)*DEG2RAD)-U(K5,K)*SIN(XC(K5)*DEG2RAD)
       VK5 = -V(K5,K)*SIN(XC(K5)*DEG2RAD)+U(K5,K)*COS(XC(K5)*DEG2RAD)
       UK6 = -V(K6,K)*COS(XC(K6)*DEG2RAD)-U(K6,K)*SIN(XC(K6)*DEG2RAD)
       VK6 = -V(K6,K)*SIN(XC(K6)*DEG2RAD)+U(K6,K)*COS(XC(K6)*DEG2RAD)

       COFA1=A1U_XY(IA,1)*UIA+A1U_XY(IA,2)*UK1   &
            +A1U_XY(IA,3)*UK2+A1U_XY(IA,4)*UK3
       COFA2=A2U_XY(IA,1)*UIA+A2U_XY(IA,2)*UK1   &
            +A2U_XY(IA,3)*UK2+A2U_XY(IA,4)*UK3
       COFA5=A1U_XY(IA,1)*VIA+A1U_XY(IA,2)*VK1   &
            +A1U_XY(IA,3)*VK2+A1U_XY(IA,4)*VK3
       COFA6=A2U_XY(IA,1)*VIA+A2U_XY(IA,2)*VK1   &
            +A2U_XY(IA,3)*VK2+A2U_XY(IA,4)*VK3

       UIJ1(II)=COFA1*XIJA+COFA2*YIJA
       VIJ1(II)=COFA5*XIJA+COFA6*YIJA
       UALFA_TMP=ABS(UIA-UIB)/ABS(UIJ1(II)+EPSILON(EPS))
       VALFA_TMP=ABS(VIA-VIB)/ABS(VIJ1(II)+EPSILON(EPS))
       IF(UALFA_TMP > 1)UALFA_TMP = 1.0_SP
       IF(VALFA_TMP > 1)VALFA_TMP = 1.0_SP
       UALFA(IA)=MIN(UALFA(IA),UALFA_TMP)
       VALFA(IA)=MIN(VALFA(IA),VALFA_TMP)

       COFA3=A1U_XY(IB,1)*UIB+A1U_XY(IB,2)*UK4   &
            +A1U_XY(IB,3)*UK5+A1U_XY(IB,4)*UK6
       COFA4=A2U_XY(IB,1)*UIB+A2U_XY(IB,2)*UK4   &
            +A2U_XY(IB,3)*UK5+A2U_XY(IB,4)*UK6
       COFA7=A1U_XY(IB,1)*VIB+A1U_XY(IB,2)*VK4   &
            +A1U_XY(IB,3)*VK5+A1U_XY(IB,4)*VK6
       COFA8=A2U_XY(IB,1)*VIB+A2U_XY(IB,2)*VK4   &
            +A2U_XY(IB,3)*VK5+A2U_XY(IB,4)*VK6

       UIJ2(II)=COFA3*XIJB+COFA4*YIJB
       VIJ2(II)=COFA7*XIJB+COFA8*YIJB
       UALFA_TMP=ABS(UIA-UIB)/ABS(UIJ2(II)+EPSILON(EPS))
       VALFA_TMP=ABS(VIA-VIB)/ABS(VIJ2(II)+EPSILON(EPS))
       IF(UALFA_TMP > 1)UALFA_TMP = 1.0_SP
       IF(VALFA_TMP > 1)VALFA_TMP = 1.0_SP
       UALFA(IB)=MIN(UALFA(IB),UALFA_TMP)
       VALFA(IB)=MIN(VALFA(IB),VALFA_TMP)

!-------ADD THE VISCOUS TERM & ADVECTION TERM---------------------------------!
!
       VISCOF1=ART(IA)*SQRT(COFA1**2+COFA6**2+0.5_SP*(COFA2+COFA5)**2)
       VISCOF2=ART(IB)*SQRT(COFA3**2+COFA8**2+0.5_SP*(COFA4+COFA7)**2)

!       VISCOF=FACT*0.5_SP*HORCON*(VISCOF1+VISCOF2)/HPRNU + FM1*HORCON
!       VISCOF=FACT*0.5_SP*HORCON*(VISCOF1+VISCOF2)/HPRNU + FM1*HORCON/HPRNU
       ! David moved HPRNU and added VHC
       VISCOF=(FACT*0.5_SP*(VISCOF1*CC_HVC(IA)+VISCOF2*CC_HVC(IB)) + FM1*0.5_SP*(CC_HVC(IA)+CC_HVC(IB)))/HPRNU

       VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD)!&
!                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))
       VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD)!&
!                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))

       VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD)!&
!                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))
       VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD)!&
!                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))

       DLTXC_TMP = VX2_TMP-VX1_TMP
       DLTYC_TMP = VY2_TMP-VY1_TMP
       
       TXXIJ=(COFA1+COFA3)*VISCOF
       TYYIJ=(COFA6+COFA8)*VISCOF
       TXYIJ=0.5_SP*(COFA2+COFA4+COFA5+COFA7)*VISCOF
       FXX(II)=DIJ*(TXXIJ*DLTYC_TMP-TXYIJ*DLTXC_TMP)
       FYY(II)=DIJ*(TXYIJ*DLTYC_TMP-TYYIJ*DLTXC_TMP)
     ENDDO

     DO II=1, NPE
       I=NPEDGE_LST(II)
       IA=IEC(I,1)
       IB=IEC(I,2)
       J1=IENODE(I,1)
       J2=IENODE(I,2)

# if !defined (SEMI_IMPLICIT)

       ELIJ=0.5_SP*(EGF(J1)+EGF(J2))
!      ggao 0103-2007   consider the sea surface pressure change
#      if defined (AIR_PRESSURE)
       ELIJ=ELIJ-0.5_SP*(EGF_AIR(J1)+EGF_AIR(J2))
#      endif
!      change end

#      if defined (EQUI_TIDE)
       ELIJ=ELIJ-0.5_SP*(EGF_EQI(J1)+EGF_EQI(J2))
#      endif
#      if defined (ATMO_TIDE)
       ELIJ=ELIJ-0.5_SP*(EGF_ATMO(J1)+EGF_ATMO(J2))
#      endif       

# else
 
       IF(STG<K_STG) THEN
         ELIJ=0.5_SP*(ET(J1)+ET(J2))
       ELSE 
         ELIJ=(1.0_SP-CETA)*0.5_SP*(ET(J1)+ET(J2))
       ENDIF
#      if defined (AIR_PRESSURE)
       ELIJ=ELIJ-( (1.0_SP-CETA)*0.5_SP*(EL_AIR(J1)+EL_AIR(J2))+CETA*0.5_SP*(ELF_AIR(J1)+ELF_AIR(J2)) )
#      endif
#      if defined (EQUI_TIDE)
       ELIJ=ELIJ-( (1.0_SP-CETA)*0.5_SP*(EL_EQI(J1)+EL_EQI(J2))+CETA*0.5_SP*(ELF_EQI(J1)+ELF_EQI(J2)) )
#      endif
#      if defined (ATMO_TIDE)
       ELIJ=ELIJ-( (1.0_SP-CETA)*0.5_SP*(EL_ATMO(J1)+EL_ATMO(J2))+CETA*0.5_SP*(ELF_ATMO(J1)+ELF_ATMO(J2)) )
#      endif

# endif

       DIJ=0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))
       UIA = -V(IA,K)*COS(XC(IA)*DEG2RAD)-U(IA,K)*SIN(XC(IA)*DEG2RAD)
       VIA = -V(IA,K)*SIN(XC(IA)*DEG2RAD)+U(IA,K)*COS(XC(IA)*DEG2RAD)
       UIB = -V(IB,K)*COS(XC(IB)*DEG2RAD)-U(IB,K)*SIN(XC(IB)*DEG2RAD)
       VIB = -V(IB,K)*SIN(XC(IB)*DEG2RAD)+U(IB,K)*COS(XC(IB)*DEG2RAD)

       UIJ1(II)=UIA+UALFA(IA)*UIJ1(II)
       VIJ1(II)=VIA+VALFA(IA)*VIJ1(II)
       UIJ2(II)=UIB+UALFA(IB)*UIJ2(II)
       VIJ2(II)=VIB+VALFA(IB)*VIJ2(II)

#      if defined (LIMITED_1)
       IF(UIJ1(II) > MAX(UIA,UIB) .OR. UIJ1(II) < MIN(UIA,UIB) .OR.  &
          UIJ2(II) > MAX(UIA,UIB) .OR. UIJ2(II) < MIN(UIA,UIB))THEN
         UIJ1(II)=UIA
         UIJ2(II)=UIB
       END IF

       IF(VIJ1(II) > MAX(VIA,VIB) .OR. VIJ1(II) < MIN(VIA,VIB) .OR.  &
          VIJ2(II) > MAX(VIA,VIB) .OR. VIJ2(II) < MIN(VIA,VIB))THEN
         VIJ1(II)=VIA
         VIJ2(II)=VIB
       END IF
#      endif       

       UIJ1_TMP = UIJ1(II)
       VIJ1_TMP = VIJ1(II)
       UIJ2_TMP = UIJ2(II)
       VIJ2_TMP = VIJ2(II)

       VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD)!&
!                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))
       VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD)!&
!                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))

       VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD)!&
!                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))
       VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD)!&
!                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))

       DLTXC_TMP = VX2_TMP-VX1_TMP
       DLTYC_TMP = VY2_TMP-VY1_TMP

       UIJ_TMP=0.5_SP*(UIJ1(II)+UIJ2(II))
       VIJ_TMP=0.5_SP*(VIJ1(II)+VIJ2(II))
       EXFLUX_TMP = DIJ*(-UIJ_TMP*DLTYC_TMP + VIJ_TMP*DLTXC_TMP)

       !!CALCULATE BOUNDARY FLUX AUGMENTERS
       TPA = FLOAT(1-ISBC(I))*EPOR(IA)
       TPB = FLOAT(1-ISBC(I))*EPOR(IB)

       !!ACCUMULATE ADVECTIVE + DIFFUSIVE + BAROTROPIC PRESSURE GRADIENT TERMS
       XADV_TMP=EXFLUX_TMP*&
                ((1.0_SP-SIGN(1.0_SP,EXFLUX_TMP))*UIJ2_TMP                 &
                +(1.0_SP+SIGN(1.0_SP,EXFLUX_TMP))*UIJ1_TMP)*0.5_SP
       YADV_TMP=EXFLUX_TMP* &
                ((1.0_SP-SIGN(1.0_SP,EXFLUX_TMP))*VIJ2_TMP                 &
                +(1.0_SP+SIGN(1.0_SP,EXFLUX_TMP))*VIJ1_TMP)*0.5_SP
       IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
#        if !defined (SEMI_IMPLICIT)
         XFLUX(IA,K)=XFLUX(IA,K)+XADV_TMP*TPA+(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX(IA,K)=YFLUX(IA,K)+YADV_TMP*TPA+(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IA)
         XFLUX(IB,K)=XFLUX(IB,K)-XADV_TMP*TPB-(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX(IB,K)=YFLUX(IB,K)-YADV_TMP*TPB-(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IB)
#        else
         XFLUX3_NP(IA,K)=XFLUX3_NP(IA,K)+XADV_TMP*TPA+(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX3_NP(IA,K)=YFLUX3_NP(IA,K)+YADV_TMP*TPA+(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IA)
         XFLUX3_NP(IB,K)=XFLUX3_NP(IB,K)-XADV_TMP*TPB-(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX3_NP(IB,K)=YFLUX3_NP(IB,K)-YADV_TMP*TPB-(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IB)
#        endif
       ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
#        if !defined (SEMI_IMPLICIT)
         XFLUX(IA,K)=XFLUX(IA,K)+XADV_TMP*TPA+(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX(IA,K)=YFLUX(IA,K)+YADV_TMP*TPA+(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IA)
#        else
         XFLUX3_NP(IA,K)=XFLUX3_NP(IA,K)+XADV_TMP*TPA+(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX3_NP(IA,K)=YFLUX3_NP(IA,K)+YADV_TMP*TPA+(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IA)
#        endif
       ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN 
#        if !defined (SEMI_IMPLICIT)
         XFLUX(IB,K)=XFLUX(IB,K)-XADV_TMP*TPB-(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX(IB,K)=YFLUX(IB,K)-YADV_TMP*TPB-(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IB)
#        else
         XFLUX3_NP(IB,K)=XFLUX3_NP(IB,K)-XADV_TMP*TPB-(FXX(II)+3.0_SP*FXX(II)*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX3_NP(IB,K)=YFLUX3_NP(IB,K)-YADV_TMP*TPB-(FYY(II)+3.0_SP*FYY(II)*FLOAT(ISBC(I)))*EPOR(IB)
#        endif
       END IF

       VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
                 * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
       VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
                 * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))

       VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
                 * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
       VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
                 * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))

       DLTXC_TMP = VX2_TMP-VX1_TMP
       DLTYC_TMP = VY2_TMP-VY1_TMP

       IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
         PSTX_TM(IA,K)=PSTX_TM(IA,K)-GRAV_E(IA)*DT1(IA)*DZ1(IA,K)*ELIJ*DLTYC_TMP/  &
                      (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
         PSTY_TM(IA,K)=PSTY_TM(IA,K)+GRAV_E(IA)*DT1(IA)*DZ1(IA,K)*ELIJ*DLTXC_TMP/  &
                      (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
         PSTX_TM(IB,K)=PSTX_TM(IB,K)+GRAV_E(IB)*DT1(IB)*DZ1(IB,K)*ELIJ*DLTYC_TMP/  &
                      (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
         PSTY_TM(IB,K)=PSTY_TM(IB,K)-GRAV_E(IB)*DT1(IB)*DZ1(IB,K)*ELIJ*DLTXC_TMP/  &
                      (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
!         PSTX_TM(IA,K)=PSTX_TM(IA,K)-GRAV_E(IA)*H1(IA)*DZ1(IA,K)*ELIJ*DLTYC_TMP/  &
!                      (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
!         PSTY_TM(IA,K)=PSTY_TM(IA,K)+GRAV_E(IA)*H1(IA)*DZ1(IA,K)*ELIJ*DLTXC_TMP/  &
!                      (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
!         PSTX_TM(IB,K)=PSTX_TM(IB,K)+GRAV_E(IB)*H1(IB)*DZ1(IB,K)*ELIJ*DLTYC_TMP/  &
!                      (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
!         PSTY_TM(IB,K)=PSTY_TM(IB,K)-GRAV_E(IB)*H1(IB)*DZ1(IB,K)*ELIJ*DLTXC_TMP/  &
!                      (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
       ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
         PSTX_TM(IA,K)=PSTX_TM(IA,K)-GRAV_E(IA)*DT1(IA)*DZ1(IA,K)*ELIJ*DLTYC_TMP/   &
                      (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
         PSTY_TM(IA,K)=PSTY_TM(IA,K)+GRAV_E(IA)*DT1(IA)*DZ1(IA,K)*ELIJ*DLTXC_TMP/   &
                      (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
!         PSTX_TM(IA,K)=PSTX_TM(IA,K)-GRAV_E(IA)*H1(IA)*DZ1(IA,K)*ELIJ*DLTYC_TMP/   &
!                      (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
!         PSTY_TM(IA,K)=PSTY_TM(IA,K)+GRAV_E(IA)*H1(IA)*DZ1(IA,K)*ELIJ*DLTXC_TMP/   &
!                      (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
       ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN       
         PSTX_TM(IB,K)=PSTX_TM(IB,K)+GRAV_E(IB)*DT1(IB)*DZ1(IB,K)*ELIJ*DLTYC_TMP/  &
                      (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
         PSTY_TM(IB,K)=PSTY_TM(IB,K)-GRAV_E(IB)*DT1(IB)*DZ1(IB,K)*ELIJ*DLTXC_TMP/  &
                      (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
!         PSTX_TM(IB,K)=PSTX_TM(IB,K)+GRAV_E(IB)*H1(IB)*DZ1(IB,K)*ELIJ*DLTYC_TMP/  &
!                      (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
!         PSTY_TM(IB,K)=PSTY_TM(IB,K)-GRAV_E(IB)*H1(IB)*DZ1(IB,K)*ELIJ*DLTXC_TMP/  &
!                      (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
       END IF
     END DO
   END DO

   DEALLOCATE(UIJ1,VIJ1,UIJ2,VIJ2,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for UIJ1,VIJ1,UIJ2 and VIJ2", &
   &                  "in subroutine ADV_UV_EDGE_XY.")
   DEALLOCATE(UALFA,VALFA,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for UALFA,VALFA", &
   &                  "in subroutine ADV_UV_EDGE_XY.")
   DEALLOCATE(FXX,FYY,STAT=ERROR)
   IF(ERROR /= 0) &
   & CALL FATAL_ERROR("Unexpected deallocation error for FXX,FYY", &
   &                  "in subroutine ADV_UV_EDGE_XY.")

#  else

   DO II=1,NPE
     I=NPEDGE_LST(II)
     IA=IEC(I,1)
     IB=IEC(I,2)
     J1=IENODE(I,1)
     J2=IENODE(I,2)
!     DIJ=0.5_SP*(DT(J1)+DT(J2))

# if !defined (SEMI_IMPLICIT)

     ELIJ=0.5_SP*(EGF(J1)+EGF(J2))
!   ggao 0103-2007   consider the sea surface pressure change
# if defined (AIR_PRESSURE)
     ELIJ=ELIJ-0.5_SP*(EGF_AIR(J1)+EGF_AIR(J2))
# endif
!     change end

#  if defined (EQUI_TIDE)
     ELIJ=ELIJ-0.5_SP*(EGF_EQI(J1)+EGF_EQI(J2))
#  endif
#  if defined (ATMO_TIDE)
     ELIJ=ELIJ-0.5_SP*(EGF_ATMO(J1)+EGF_ATMO(J2))
#    endif       

# else
 
     IF(STG<K_STG) THEN
       ELIJ=0.5_SP*(ET(J1)+ET(J2))
     ELSE
       ELIJ=(1.0_SP-CETA)*0.5_SP*(ET(J1)+ET(J2))
     ENDIF
#    if defined (AIR_PRESSURE)
     ELIJ=ELIJ-( (1.0_SP-CETA)*0.5_SP*(EL_AIR(J1)+EL_AIR(J2))+CETA*0.5_SP*(ELF_AIR(J1)+ELF_AIR(J2)) )
#    endif
#    if defined (EQUI_TIDE)
     ELIJ=ELIJ-( (1.0_SP-CETA)*0.5_SP*(EL_EQI(J1)+EL_EQI(J2))+CETA*0.5_SP*(ELF_EQI(J1)+ELF_EQI(J2)) )
#    endif
#    if defined (ATMO_TIDE)
     ELIJ=ELIJ-( (1.0_SP-CETA)*0.5_SP*(EL_ATMO(J1)+EL_ATMO(J2))+CETA*0.5_SP*(ELF_ATMO(J1)+ELF_ATMO(J2)) )
#    endif

# endif

     K1=NBE(IA,1)
     K2=NBE(IA,2)
     K3=NBE(IA,3)
     K4=NBE(IB,1)
     K5=NBE(IB,2)
     K6=NBE(IB,3)

     XIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * COS(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
     YIJC_TMP = REARTH * COS(YIJC(I)*DEG2RAD) * SIN(XIJC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YIJC(I)*DEG2RAD))
		  
     XCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * COS(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
     YCIA_TMP = REARTH * COS(YC(IA)*DEG2RAD) * SIN(XC(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IA)*DEG2RAD))
     XCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * COS(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
     YCIB_TMP = REARTH * COS(YC(IB)*DEG2RAD) * SIN(XC(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+SIN(YC(IB)*DEG2RAD))
		  
     XIJA = XIJC_TMP-XCIA_TMP
     YIJA = YIJC_TMP-YCIA_TMP
     XIJB = XIJC_TMP-XCIB_TMP
     YIJB = YIJC_TMP-YCIB_TMP

     DO K=1,KBM1

       DIJ=0.5_SP*(DT(J1)*DZ(J1,K)+DT(J2)*DZ(J2,K))

       UIA = -V(IA,K)*COS(XC(IA)*DEG2RAD)-U(IA,K)*SIN(XC(IA)*DEG2RAD)
       VIA = -V(IA,K)*SIN(XC(IA)*DEG2RAD)+U(IA,K)*COS(XC(IA)*DEG2RAD)
       UIB = -V(IB,K)*COS(XC(IB)*DEG2RAD)-U(IB,K)*SIN(XC(IB)*DEG2RAD)
       VIB = -V(IB,K)*SIN(XC(IB)*DEG2RAD)+U(IB,K)*COS(XC(IB)*DEG2RAD)
       UK1 = -V(K1,K)*COS(XC(K1)*DEG2RAD)-U(K1,K)*SIN(XC(K1)*DEG2RAD)
       VK1 = -V(K1,K)*SIN(XC(K1)*DEG2RAD)+U(K1,K)*COS(XC(K1)*DEG2RAD)
       UK2 = -V(K2,K)*COS(XC(K2)*DEG2RAD)-U(K2,K)*SIN(XC(K2)*DEG2RAD)
       VK2 = -V(K2,K)*SIN(XC(K2)*DEG2RAD)+U(K2,K)*COS(XC(K2)*DEG2RAD)
       UK3 = -V(K3,K)*COS(XC(K3)*DEG2RAD)-U(K3,K)*SIN(XC(K3)*DEG2RAD)
       VK3 = -V(K3,K)*SIN(XC(K3)*DEG2RAD)+U(K3,K)*COS(XC(K3)*DEG2RAD)
       UK4 = -V(K4,K)*COS(XC(K4)*DEG2RAD)-U(K4,K)*SIN(XC(K4)*DEG2RAD)
       VK4 = -V(K4,K)*SIN(XC(K4)*DEG2RAD)+U(K4,K)*COS(XC(K4)*DEG2RAD)
       UK5 = -V(K5,K)*COS(XC(K5)*DEG2RAD)-U(K5,K)*SIN(XC(K5)*DEG2RAD)
       VK5 = -V(K5,K)*SIN(XC(K5)*DEG2RAD)+U(K5,K)*COS(XC(K5)*DEG2RAD)
       UK6 = -V(K6,K)*COS(XC(K6)*DEG2RAD)-U(K6,K)*SIN(XC(K6)*DEG2RAD)
       VK6 = -V(K6,K)*SIN(XC(K6)*DEG2RAD)+U(K6,K)*COS(XC(K6)*DEG2RAD)

       COFA1=A1U_XY(IA,1)*UIA+A1U_XY(IA,2)*UK1   &
            +A1U_XY(IA,3)*UK2+A1U_XY(IA,4)*UK3
       COFA2=A2U_XY(IA,1)*UIA+A2U_XY(IA,2)*UK1   &
            +A2U_XY(IA,3)*UK2+A2U_XY(IA,4)*UK3
       COFA5=A1U_XY(IA,1)*VIA+A1U_XY(IA,2)*VK1   &
            +A1U_XY(IA,3)*VK2+A1U_XY(IA,4)*VK3
       COFA6=A2U_XY(IA,1)*VIA+A2U_XY(IA,2)*VK1   &
            +A2U_XY(IA,3)*VK2+A2U_XY(IA,4)*VK3

       UIJ1=UIA+COFA1*XIJA+COFA2*YIJA
       VIJ1=VIA+COFA5*XIJA+COFA6*YIJA

       COFA3=A1U_XY(IB,1)*UIB+A1U_XY(IB,2)*UK4   &
            +A1U_XY(IB,3)*UK5+A1U_XY(IB,4)*UK6
       COFA4=A2U_XY(IB,1)*UIB+A2U_XY(IB,2)*UK4   &
            +A2U_XY(IB,3)*UK5+A2U_XY(IB,4)*UK6
       COFA7=A1U_XY(IB,1)*VIB+A1U_XY(IB,2)*VK4   &
            +A1U_XY(IB,3)*VK5+A1U_XY(IB,4)*VK6
       COFA8=A2U_XY(IB,1)*VIB+A2U_XY(IB,2)*VK4   &
            +A2U_XY(IB,3)*VK5+A2U_XY(IB,4)*VK6

       UIJ2=UIB+COFA3*XIJB+COFA4*YIJB
       VIJ2=VIB+COFA7*XIJB+COFA8*YIJB

!-------ADD THE VISCOUS TERM & ADVECTION TERM---------------------------------!
!
       VISCOF1=ART(IA)*SQRT(COFA1**2+COFA6**2+0.5_SP*(COFA2+COFA5)**2)
       VISCOF2=ART(IB)*SQRT(COFA3**2+COFA8**2+0.5_SP*(COFA4+COFA7)**2)

       !VISCOF=FACT*0.5_SP*HORCON*(VISCOF1+VISCOF2)/HPRNU + FM1*HORCON
     
       ! David moved HPRNU and added VHC
       VISCOF=(FACT*0.5_SP*(VISCOF1*CC_HVC(IA)+VISCOF2*CC_HVC(IB)) + FM1*0.5_SP*(CC_HVC(IA)+CC_HVC(IB)))/HPRNU

       VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD)       &
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))
       VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,1))*DEG2RAD))

       VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))
       VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD)&
                  * 2._SP /(1._SP+SIN(VY(IENODE(I,2))*DEG2RAD))

       DLTXC_TMP = VX2_TMP-VX1_TMP
       DLTYC_TMP = VY2_TMP-VY1_TMP
       
       TXXIJ=(COFA1+COFA3)*VISCOF
       TYYIJ=(COFA6+COFA8)*VISCOF
       TXYIJ=0.5_SP*(COFA2+COFA4+COFA5+COFA7)*VISCOF
       FXX=DIJ*(TXXIJ*DLTYC_TMP-TXYIJ*DLTXC_TMP)
       FYY=DIJ*(TXYIJ*DLTYC_TMP-TYYIJ*DLTXC_TMP)


       UIJ1_TMP = UIJ1
       VIJ1_TMP = VIJ1
       UIJ2_TMP = UIJ2
       VIJ2_TMP = VIJ2

       UIJ_TMP=0.5_SP*(UIJ1+UIJ2)
       VIJ_TMP=0.5_SP*(VIJ1+VIJ2)
       EXFLUX_TMP = DIJ*(-UIJ_TMP*DLTYC_TMP + VIJ_TMP*DLTXC_TMP)

       !!CALCULATE BOUNDARY FLUX AUGMENTERS
       TPA = FLOAT(1-ISBC(I))*EPOR(IA)
       TPB = FLOAT(1-ISBC(I))*EPOR(IB)

       !!ACCUMULATE ADVECTIVE + DIFFUSIVE + BAROTROPIC PRESSURE GRADIENT TERMS
       XADV_TMP=EXFLUX_TMP*&
                ((1.0_SP-SIGN(1.0_SP,EXFLUX_TMP))*UIJ2_TMP                 &
                +(1.0_SP+SIGN(1.0_SP,EXFLUX_TMP))*UIJ1_TMP)*0.5_SP
       YADV_TMP=EXFLUX_TMP* &
                ((1.0_SP-SIGN(1.0_SP,EXFLUX_TMP))*VIJ2_TMP                 &
 	        +(1.0_SP+SIGN(1.0_SP,EXFLUX_TMP))*VIJ1_TMP)*0.5_SP
       IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
#        if !defined (SEMI_IMPLICIT)
         XFLUX(IA,K)=XFLUX(IA,K)+XADV_TMP*TPA+(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX(IA,K)=YFLUX(IA,K)+YADV_TMP*TPA+(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IA)
         XFLUX(IB,K)=XFLUX(IB,K)-XADV_TMP*TPB-(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX(IB,K)=YFLUX(IB,K)-YADV_TMP*TPB-(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IB)
#        else
         XFLUX3_NP(IA,K)=XFLUX3_NP(IA,K)+XADV_TMP*TPA+(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX3_NP(IA,K)=YFLUX3_NP(IA,K)+YADV_TMP*TPA+(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IA)
         XFLUX3_NP(IB,K)=XFLUX3_NP(IB,K)-XADV_TMP*TPB-(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX3_NP(IB,K)=YFLUX3_NP(IB,K)-YADV_TMP*TPB-(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IB)
#        endif
       ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
#        if !defined (SEMI_IMPLICIT)
         XFLUX(IA,K)=XFLUX(IA,K)+XADV_TMP*TPA+(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX(IA,K)=YFLUX(IA,K)+YADV_TMP*TPA+(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IA)
#        else
         XFLUX3_NP(IA,K)=XFLUX3_NP(IA,K)+XADV_TMP*TPA+(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IA)
         YFLUX3_NP(IA,K)=YFLUX3_NP(IA,K)+YADV_TMP*TPA+(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IA)
#        endif
       ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN 
#        if !defined (SEMI_IMPLICIT)
         XFLUX(IB,K)=XFLUX(IB,K)-XADV_TMP*TPB-(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX(IB,K)=YFLUX(IB,K)-YADV_TMP*TPB-(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IB)
#        else
         XFLUX3_NP(IB,K)=XFLUX3_NP(IB,K)-XADV_TMP*TPB-(FXX+3.0_SP*FXX*FLOAT(ISBC(I)))*EPOR(IB)
         YFLUX3_NP(IB,K)=YFLUX3_NP(IB,K)-YADV_TMP*TPB-(FYY+3.0_SP*FYY*FLOAT(ISBC(I)))*EPOR(IB)
#        endif
       END IF

       VX1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * COS(VX(IENODE(I,1))*DEG2RAD) &
                 * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))
       VY1_TMP = REARTH * COS(VY(IENODE(I,1))*DEG2RAD) * SIN(VX(IENODE(I,1))*DEG2RAD) &
                 * 2._SP /(1._SP+sin(VY(IENODE(I,1))*DEG2RAD))

       VX2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * COS(VX(IENODE(I,2))*DEG2RAD) &
                 * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))
       VY2_TMP = REARTH * COS(VY(IENODE(I,2))*DEG2RAD) * SIN(VX(IENODE(I,2))*DEG2RAD) &
                 * 2._SP /(1._SP+sin(VY(IENODE(I,2))*DEG2RAD))

       DLTXC_TMP = VX2_TMP-VX1_TMP
       DLTYC_TMP = VY2_TMP-VY1_TMP

       IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) == 1)THEN
         PSTX_TM(IA,K)=PSTX_TM(IA,K)-GRAV_E(IA)*DT1(IA)*DZ1(IA,K)*ELIJ*DLTYC_TMP/  &
	              (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
         PSTY_TM(IA,K)=PSTY_TM(IA,K)+GRAV_E(IA)*DT1(IA)*DZ1(IA,K)*ELIJ*DLTXC_TMP/  &
	              (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
         PSTX_TM(IB,K)=PSTX_TM(IB,K)+GRAV_E(IB)*DT1(IB)*DZ1(IB,K)*ELIJ*DLTYC_TMP/  &
	              (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
         PSTY_TM(IB,K)=PSTY_TM(IB,K)-GRAV_E(IB)*DT1(IB)*DZ1(IB,K)*ELIJ*DLTXC_TMP/  &
	              (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
       ELSE IF(CELL_NORTHAREA(IA) == 1 .AND. CELL_NORTHAREA(IB) /= 1)THEN
         PSTX_TM(IA,K)=PSTX_TM(IA,K)-GRAV_E(IA)*DT1(IA)*DZ1(IA,K)*ELIJ*DLTYC_TMP/   &
	              (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
         PSTY_TM(IA,K)=PSTY_TM(IA,K)+GRAV_E(IA)*DT1(IA)*DZ1(IA,K)*ELIJ*DLTXC_TMP/   &
	              (2._SP /(1._SP+sin(YC(IA)*DEG2RAD)))
       ELSE IF(CELL_NORTHAREA(IB) == 1 .AND. CELL_NORTHAREA(IA) /= 1)THEN       
         PSTX_TM(IB,K)=PSTX_TM(IB,K)+GRAV_E(IB)*DT1(IB)*DZ1(IB,K)*ELIJ*DLTYC_TMP/  &
	              (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
         PSTY_TM(IB,K)=PSTY_TM(IB,K)-GRAV_E(IB)*DT1(IB)*DZ1(IB,K)*ELIJ*DLTXC_TMP/  &
	              (2._SP /(1._SP+sin(YC(IB)*DEG2RAD)))
       END IF
     END DO
   END DO

#  endif

   DO K = 1,KBM1
     DO II=1,NP
       I=NP_LST(II)
#      if !defined (SEMI_IMPLICIT)
       XFLUX(I,K)=XFLUX(I,K)+PSTX_TM(I,K)
       YFLUX(I,K)=YFLUX(I,K)+PSTY_TM(I,K)
#      else
       XFLUX3_NP(I,K)=XFLUX3_NP(I,K)+PSTX_TM(I,K)
       YFLUX3_NP(I,K)=YFLUX3_NP(I,K)+PSTY_TM(I,K)
#      endif
     END DO
   END DO

!
!-------ADD VERTICAL CONVECTIVE FLUX, CORIOLIS TERM AND BAROCLINIC PG TERM----!
!
   DO II=1,NP
     I=NP_LST(II)
     IF(CELL_NORTHAREA(I) == 1)THEN
       DO K=1,KBM1
         IF(K == 1) THEN
           UIJ1_TMP = -V(I,K)*COS(XC(I)*DEG2RAD)-U(I,K)*SIN(XC(I)*DEG2RAD)
           VIJ1_TMP = -V(I,K)*SIN(XC(I)*DEG2RAD)+U(I,K)*COS(XC(I)*DEG2RAD)
           UIJ2_TMP = -V(I,K+1)*COS(XC(I)*DEG2RAD)-U(I,K+1)*SIN(XC(I)*DEG2RAD)
           VIJ2_TMP = -V(I,K+1)*SIN(XC(I)*DEG2RAD)+U(I,K+1)*COS(XC(I)*DEG2RAD)
           XFLUXV=-W(I,K+1)*(UIJ1_TMP*DZ1(I,K+1)+UIJ2_TMP*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K+1))
           YFLUXV=-W(I,K+1)*(VIJ1_TMP*DZ1(I,K+1)+VIJ2_TMP*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K+1))
         ELSE IF(K == KBM1) THEN
           UIJ1_TMP = -V(I,K)*COS(XC(I)*DEG2RAD)-U(I,K)*SIN(XC(I)*DEG2RAD)
           VIJ1_TMP = -V(I,K)*SIN(XC(I)*DEG2RAD)+U(I,K)*COS(XC(I)*DEG2RAD)
           UIJ2_TMP = -V(I,K-1)*COS(XC(I)*DEG2RAD)-U(I,K-1)*SIN(XC(I)*DEG2RAD)
           VIJ2_TMP = -V(I,K-1)*SIN(XC(I)*DEG2RAD)+U(I,K-1)*COS(XC(I)*DEG2RAD)
           XFLUXV= W(I,K)*(UIJ1_TMP*DZ1(I,K-1)+UIJ2_TMP*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K-1))
           YFLUXV= W(I,K)*(VIJ1_TMP*DZ1(I,K-1)+VIJ2_TMP*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K-1))
         ELSE
           UIJ1_TMP = -V(I,K)*COS(XC(I)*DEG2RAD)-U(I,K)*SIN(XC(I)*DEG2RAD)
           VIJ1_TMP = -V(I,K)*SIN(XC(I)*DEG2RAD)+U(I,K)*COS(XC(I)*DEG2RAD)
           UIJ2_TMP = -V(I,K-1)*COS(XC(I)*DEG2RAD)-U(I,K-1)*SIN(XC(I)*DEG2RAD)
           VIJ2_TMP = -V(I,K-1)*SIN(XC(I)*DEG2RAD)+U(I,K-1)*COS(XC(I)*DEG2RAD)
           UIJ3_TMP = -V(I,K+1)*COS(XC(I)*DEG2RAD)-U(I,K+1)*SIN(XC(I)*DEG2RAD)
           VIJ3_TMP = -V(I,K+1)*SIN(XC(I)*DEG2RAD)+U(I,K+1)*COS(XC(I)*DEG2RAD)
           XFLUXV= W(I,K)*(UIJ1_TMP*DZ1(I,K-1)+UIJ2_TMP*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K-1))-  &
                   W(I,K+1)*(UIJ1_TMP*DZ1(I,K+1)+UIJ3_TMP*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K+1))
           YFLUXV= W(I,K)*(VIJ1_TMP*DZ1(I,K-1)+VIJ2_TMP*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K-1))-  &
                   W(I,K+1)*(VIJ1_TMP*DZ1(I,K+1)+VIJ3_TMP*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K+1))
         END IF
         U_TMP = -V(I,K)*COS(XC(I)*DEG2RAD)-U(I,K)*SIN(XC(I)*DEG2RAD)
         V_TMP = -V(I,K)*SIN(XC(I)*DEG2RAD)+U(I,K)*COS(XC(I)*DEG2RAD)
#        if !defined (SEMI_IMPLICIT)
!         XFLUX(I,K)=XFLUX(I,K)+XFLUXV/DZ(K)*ART(I)&
!                    +DRHOX(I,K)-COR(I)*V_TMP*DT1(I)*ART(I)
!         YFLUX(I,K)=YFLUX(I,K)+YFLUXV/DZ(K)*ART(I)&
!                    +DRHOY(I,K)+COR(I)*U_TMP*DT1(I)*ART(I)
         XFLUX(I,K)=XFLUX(I,K)+XFLUXV*ART(I)&
                    +DRHOX(I,K)-COR(I)*V_TMP*DT1(I)*DZ1(I,K)*ART(I)
         YFLUX(I,K)=YFLUX(I,K)+YFLUXV*ART(I)&
                    +DRHOY(I,K)+COR(I)*U_TMP*DT1(I)*DZ1(I,K)*ART(I)
#        else
!         XFLUX3_NP(I,K)=XFLUX3_NP(I,K)+XFLUXV/DZ(K)*ART(I)&
!                    +DRHOX(I,K)-COR(I)*V_TMP*DT1(I)*ART(I)
!         YFLUX3_NP(I,K)=YFLUX3_NP(I,K)+YFLUXV/DZ(K)*ART(I)&
!                    +DRHOY(I,K)+COR(I)*U_TMP*DT1(I)*ART(I)
         XFLUX3_NP(I,K)=XFLUX3_NP(I,K)+XFLUXV*ART(I)&
                    +DRHOX(I,K)-COR(I)*V_TMP*DT1(I)*DZ1(I,K)*ART(I)
         YFLUX3_NP(I,K)=YFLUX3_NP(I,K)+YFLUXV*ART(I)&
                    +DRHOY(I,K)+COR(I)*U_TMP*DT1(I)*DZ1(I,K)*ART(I)
#        endif

       END DO
     END IF
   END DO

#  if defined (SEMI_IMPLICIT)
   DO II=1, NP
     I = NP_LST(II)
     IF(CELL_NORTHAREA(I) ==1 ) THEN
       DO K=1, KBM1
         XFLUX(I,K) = YFLUX3_NP(I,K)*COS(XC(I)*DEG2RAD)-XFLUX3_NP(I,K)*SIN(XC(I)*DEG2RAD)
         YFLUX(I,K) = -( XFLUX3_NP(I,K)*COS(XC(I)*DEG2RAD)+YFLUX3_NP(I,K)*SIN(XC(I)*DEG2RAD) )         
       ENDDO
     ENDIF
   ENDDO   
#  endif

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: adv_uv_edge_xy"

   RETURN
   END SUBROUTINE ADV_UV_EDGE_XY
!==============================================================================!

!==============================================================================|
#  if !defined (SEMI_IMPLICIT)
   SUBROUTINE EXTEL_EDGE_XY(K,XFLUX)       

   USE MOD_OBCS
   IMPLICIT NONE
   REAL(SP) :: XFLUX(0:MT)
   REAL(SP) :: DIJ,UIJ,VIJ
   INTEGER  :: I,J,K,I1,IA,IB,JJ,J1,J2,II

   REAL(SP) :: UIJ_TMP,VIJ_TMP,EXFLUX_TMP
   REAL(SP) :: DLTXE_TMP,DLTYE_TMP
   REAL(SP) :: VX1_TMP,VX2_TMP,VY1_TMP,VY2_TMP

   ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN
   
   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: extel_edge_XY"
!----------INITIALIZE FLUX ARRAY ----------------------------------------------!

   DO II=1,NPCV
     I = NCEDGE_LST(II)
     IA  = NIEC(I,1)
     IB  = NIEC(I,2)
     IF(IA == NODE_NORTHPOLE)THEN
       XFLUX(IA) = 0.0_SP
     END IF
     IF(IB == NODE_NORTHPOLE)THEN  
       XFLUX(IB) = 0.0_SP
     END IF  
   END DO  

!---------ACCUMULATE FLUX BY LOOPING OVER CONTROL VOLUME HALF EDGES------------!

   DO II=1,NPCV
     I = NCEDGE_LST(II)
     I1  = NTRG(I)
     IA  = NIEC(I,1)
     IB  = NIEC(I,2)
     DIJ = D1(I1)

     UIJ = UA(I1)
     VIJ = VA(I1)

     IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
       UIJ_TMP = -VIJ*COS(XC(I1)*DEG2RAD)-UIJ*SIN(XC(I1)*DEG2RAD)
       VIJ_TMP = -VIJ*SIN(XC(I1)*DEG2RAD)+UIJ*COS(XC(I1)*DEG2RAD)
       
       VX1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD)&
                 * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
       VY1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD)&
                 * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))

       VX2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD)&
                 * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       VY2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD)&
                 * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))

       DLTXE_TMP = VX2_TMP-VX1_TMP
       DLTYE_TMP = VY2_TMP-VY1_TMP
       
       EXFLUX_TMP = DIJ*(-UIJ_TMP*DLTYE_TMP+VIJ_TMP*DLTXE_TMP)
     END IF  
     
     IF(IA == NODE_NORTHPOLE) THEN    
       XFLUX(IA) = XFLUX(IA)-EXFLUX_TMP
     ELSE IF(IB == NODE_NORTHPOLE)THEN
       XFLUX(IB) = XFLUX(IB)+EXFLUX_TMP
     END IF    

   END DO

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: extel_edge_XY"

   RETURN
   END SUBROUTINE EXTEL_EDGE_XY
#  endif
!==============================================================================|

!==============================================================================|
#  if !defined (SEMI_IMPLICIT)
   SUBROUTINE EXTUV_EDGE_XY(K)       

   IMPLICIT NONE
   INTEGER, INTENT(IN) :: K
   REAL(SP), DIMENSION(0:NT) :: RESX,RESY
   INTEGER  :: I,II

   REAL(SP) :: UARK_TMP,VARK_TMP,UAF_TMP,VAF_TMP,UA_TMP,VA_TMP

   REAL(SP) :: WUSURF2_TMP,WVSURF2_TMP,WUBOT_TMP,WVBOT_TMP
   ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN
   
   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: extuv_edge_XY"

!
!--ACCUMULATE RESIDUALS FOR EXTERNAL MODE EQUATIONS----------------------------|
!
   DO II=1,NP
     I=NP_LST(II)
     UA_TMP = -VA(I)*COS(XC(I)*DEG2RAD)-UA(I)*SIN(XC(I)*DEG2RAD)
     VA_TMP = -VA(I)*SIN(XC(I)*DEG2RAD)+UA(I)*COS(XC(I)*DEG2RAD)

     WUSURF2_TMP = -WVSURF2(I)*COS(XC(I)*DEG2RAD)-WUSURF2(I)*SIN(XC(I)*DEG2RAD)
     WVSURF2_TMP = -WVSURF2(I)*SIN(XC(I)*DEG2RAD)+WUSURF2(I)*COS(XC(I)*DEG2RAD)

     WUBOT_TMP = -WVBOT(I)*COS(XC(I)*DEG2RAD)-WUBOT(I)*SIN(XC(I)*DEG2RAD)
     WVBOT_TMP = -WVBOT(I)*SIN(XC(I)*DEG2RAD)+WUBOT(I)*COS(XC(I)*DEG2RAD)

     RESX(I) = ADX2D(I)+ADVUA(I)+DRX2D(I)+PSTX(I)                &
               -COR(I)*VA_TMP*D1(I)*ART(I)                       &
               -(WUSURF2_TMP+WUBOT_TMP)*ART(I)
     RESY(I) = ADY2D(I)+ADVVA(I)+DRY2D(I)+PSTY(I)                &
               +COR(I)*UA_TMP*D1(I)*ART(I)                       &
               -(WVSURF2_TMP+WVBOT_TMP)*ART(I)

!
!--UPDATE----------------------------------------------------------------------|
!
     UARK_TMP = -VARK(I)*COS(XC(I)*DEG2RAD)-UARK(I)*SIN(XC(I)*DEG2RAD)
     VARK_TMP = -VARK(I)*SIN(XC(I)*DEG2RAD)+UARK(I)*COS(XC(I)*DEG2RAD)

     UAF_TMP = (UARK_TMP*(H1(I)+ELRK1(I))              &
               -ALPHA_RK(K)*DTE*RESX(I)/ART(I))/(H1(I)+ELF1(I))
     VAF_TMP = (VARK_TMP*(H1(I)+ELRK1(I))              &
               -ALPHA_RK(K)*DTE*RESY(I)/ART(I))/(H1(I)+ELF1(I))
		       
     UAF(I)  = VAF_TMP*COS(XC(I)*DEG2RAD)-UAF_TMP*SIN(XC(I)*DEG2RAD)
     VAF(I)  = UAF_TMP*COS(XC(I)*DEG2RAD)+VAF_TMP*SIN(XC(I)*DEG2RAD)
     VAF(I)  = -VAF(I)			    

   END DO

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: extuv_edge_XY"

   RETURN
   END SUBROUTINE EXTUV_EDGE_XY
#  endif
!==============================================================================|

!==============================================================================|

   SUBROUTINE VERTVL_EDGE_XY(XFLUX,CETA)         

!------------------------------------------------------------------------------|
   IMPLICIT NONE 
   REAL(SP) :: XFLUX(MT,KBM1)
   REAL(SP) :: DIJ,UIJ,VIJ,UN,EXFLUX,TMP1,DIJ1,UIJ1,VIJ1
   INTEGER  :: I,K,IA,IB,I1 ,J,JJ,J1,J2,II

   REAL(SP) :: UIJ_TMP,VIJ_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,UIJ1_TMP,VIJ1_TMP
   REAL(SP) :: DLTXE_TMP,DLTYE_TMP,EXFLUX_TMP

   REAL(SP) :: CETA
!------------------------------------------------------------------------------|
   ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN
   
   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: Vertvl_edge_XY"

!----------------------INITIALIZE FLUX-----------------------------------------!

   DO K=1,KBM1
     DO II=1,NPCV
       I = NCEDGE_LST(II)
       IA  = NIEC(I,1)
       IB  = NIEC(I,2)
       IF(IA == NODE_NORTHPOLE)THEN
         XFLUX(IA,K) = 0.0_SP
       END IF
       IF(IB == NODE_NORTHPOLE)THEN  
         XFLUX(IB,K) = 0.0_SP
       END IF  
     END DO  
   END DO
!----------------------ACCUMULATE FLUX-----------------------------------------!

   DO II=1,NPCV
     I=NCEDGE_LST(II)
     I1=NTRG(I)
     IA=NIEC(I,1)
     IB=NIEC(I,2)
!     DIJ=DT1(I1)
     DO K=1,KBM1
#      if !defined (SEMI_IMPLICIT)
       DIJ=DT1(I1)*DZ1(I1,K)
       UIJ=U(I1,K)
       VIJ=V(I1,K)
#      else
       DIJ=DT1(I1)*DZ1(I1,K)
       DIJ1=D1(I1)*DZ1(I1,K)
       UIJ=U(I1,K)
       VIJ=V(I1,K)
       UIJ1=UF(I1,K)
       VIJ1=VF(I1,K)
#      endif

       IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
         UIJ_TMP = -VIJ*COS(XC(I1)*DEG2RAD)-UIJ*SIN(XC(I1)*DEG2RAD)
         VIJ_TMP = -VIJ*SIN(XC(I1)*DEG2RAD)+UIJ*COS(XC(I1)*DEG2RAD)
#        if defined (SEMI_IMPLICIT)
         UIJ1_TMP = -VIJ1*COS(XC(I1)*DEG2RAD)-UIJ1*SIN(XC(I1)*DEG2RAD)
         VIJ1_TMP = -VIJ1*SIN(XC(I1)*DEG2RAD)+UIJ1*COS(XC(I1)*DEG2RAD)
#        endif

       VX1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD)&
                 * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
       VY1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD)&
                 * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))

       VX2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD)&
                 * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       VY2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD)&
                 * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))

         DLTXE_TMP = VX2_TMP-VX1_TMP
         DLTYE_TMP = VY2_TMP-VY1_TMP
       
#        if !defined (SEMI_IMPLICIT)
         EXFLUX_TMP = DIJ*(-UIJ_TMP*DLTYE_TMP+VIJ_TMP*DLTXE_TMP)
#        else
         EXFLUX_TMP = (1.0_SP-CETA)*DIJ*(-UIJ_TMP*DLTYE_TMP+VIJ_TMP*DLTXE_TMP)+CETA*DIJ1*(-UIJ1_TMP*DLTYE_TMP+VIJ1_TMP*DLTXE_TMP)
#        endif
       END IF  
     
       IF(IA == NODE_NORTHPOLE)THEN
         XFLUX(IA,K) = XFLUX(IA,K)-EXFLUX_TMP
       ELSE IF(IB == NODE_NORTHPOLE)THEN
         XFLUX(IB,K) = XFLUX(IB,K)+EXFLUX_TMP
       END IF
     END DO
   END DO

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: Vertvl_edge_xy"

   RETURN
   END SUBROUTINE VERTVL_EDGE_XY
!==============================================================================|

!==============================================================================|
   SUBROUTINE ADV_S_XY(XFLUX,XFLUX_ADV,PSPX,PSPY,PSPXD,PSPYD,VISCOFF,K,CETA)               

!------------------------------------------------------------------------------|

   IMPLICIT NONE
   INTEGER, INTENT(IN) :: K

   REAL(SP), DIMENSION(0:MT,KB)     :: XFLUX,XFLUX_ADV
   REAL(SP), DIMENSION(M)           :: PSPX,PSPY,PSPXD,PSPYD,VISCOFF
!   REAL(SP), DIMENSION(3*(NT))      :: DTIJ 
   REAL(SP), DIMENSION(3*(NT),KBM1)      :: DTIJ 
   REAL(SP) :: XI,YI
   REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2 
   REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF   
   REAL(SP) :: FACT,FM1
   INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,JTMP,JJ,II
   REAL(SP) :: TXPI,TYPI

   REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
   REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
   REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
   REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
   REAL(SP) :: PSPX_TMP,PSPY_TMP,PSPXD_TMP,PSPYD_TMP
   REAL(SP) :: U_TMP,V_TMP
   REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2

!!  ggao edge calculation
   REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
   REAL(SP) :: S1MIN, S1MAX, S2MIN, S2MAX

   REAL(SP) :: CETA
#  if defined (SEMI_IMPLICIT)
   REAL(SP), DIMENSION(3*(NT),KBM1)  :: DTIJ1
   REAL(SP) :: UIJ_TMP1, VIJ_TMP1, UVN_TMP1
#  endif
!------------------------------------------------------------------------------!
   ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN
   
   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: ADV_S_XY(K):",K

   
  SELECT CASE(HORIZONTAL_MIXING_TYPE)
  CASE ('closure')
     FACT = 1.0_SP
     FM1  = 0.0_SP
  CASE('constant')
     FACT = 0.0_SP
     FM1  = 1.0_SP
  CASE DEFAULT
     CALL FATAL_ERROR("UNKNOW HORIZONTAL MIXING TYPE:",&
          & TRIM(HORIZONTAL_MIXING_TYPE) )
  END SELECT

!
!--Initialize Fluxes-----------------------------------------------------------!
!
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     IA = NIEC(I,1)
     IB = NIEC(I,2)
     IF(IA == NODE_NORTHPOLE)THEN
       XFLUX(IA,K) = 0.0_SP
       XFLUX_ADV(IA,K) = 0.0_SP
     ELSE IF(IB == NODE_NORTHPOLE)THEN  
       XFLUX(IB,K) = 0.0_SP
       XFLUX_ADV(IB,K) = 0.0_SP
     END IF  
   END DO  
     
!
!--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
!
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     I1=NTRG(I)
     DTIJ(I,K)=DT1(I1)*DZ1(I1,K)
#    if defined (SEMI_IMPLICIT)
     DTIJ1(I,K) = D1(I1)*DZ1(I1,K)
#    endif
   END DO

!
!--Calculate the Advection and Horizontal Diffusion Terms----------------------!
!
   I = NODE_NORTHPOLE

   IF(I==0)  RETURN

   PUPX_TMP=0.0_SP
   PUPY_TMP=0.0_SP
   PVPX_TMP=0.0_SP
   PVPY_TMP=0.0_SP

   DO J=1,NTVE(I)
     I1=NBVE(I,J)
     JTMP=NBVT(I,J)
     J1=JTMP+1-(JTMP+1)/4*3
     J2=JTMP+2-(JTMP+2)/4*3
       
     VX_TMP = REARTH * COS(VY(I)*DEG2RAD) * COS(VX(I)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(I)*DEG2RAD))
     VY_TMP = REARTH * COS(VY(I)*DEG2RAD) * SIN(VX(I)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(I)*DEG2RAD))
		     
     VX1_TMP= REARTH * COS(VY(NV(I1,J1))*DEG2RAD) * COS(VX(NV(I1,J1))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J1))*DEG2RAD))
     VY1_TMP= REARTH * COS(VY(NV(I1,J1))*DEG2RAD) * SIN(VX(NV(I1,J1))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J1))*DEG2RAD))
		     
     VX2_TMP= REARTH * COS(YC(I1)*DEG2RAD) * COS(XC(I1)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(YC(I1)*DEG2RAD))
     VY2_TMP= REARTH * COS(YC(I1)*DEG2RAD) * SIN(XC(I1)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(YC(I1)*DEG2RAD))
		     
     VX3_TMP= REARTH * COS(VY(NV(I1,J2))*DEG2RAD) * COS(VX(NV(I1,J2))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J2))*DEG2RAD))
     VY3_TMP= REARTH * COS(VY(NV(I1,J2))*DEG2RAD) * SIN(VX(NV(I1,J2))*DEG2RAD) &
                    * 2._SP /(1._SP+SIN(VY(NV(I1,J2))*DEG2RAD))
		     
     X11=0.5_SP*(VX_TMP+VX1_TMP)
     Y11=0.5_SP*(VY_TMP+VY1_TMP)
     X22=VX2_TMP
     Y22=VX2_TMP
     X33=0.5_SP*(VX_TMP+VX3_TMP)
     Y33=0.5_SP*(VY_TMP+VY3_TMP)
     
     U_TMP = -V(I1,K)*COS(XC(I1)*DEG2RAD)-U(I1,K)*SIN(XC(I1)*DEG2RAD)
     V_TMP = -V(I1,K)*SIN(XC(I1)*DEG2RAD)+U(I1,K)*COS(XC(I1)*DEG2RAD)

     PUPX_TMP=PUPX_TMP+U_TMP*(Y11-Y33)
     PUPY_TMP=PUPY_TMP+U_TMP*(X33-X11)
     PVPX_TMP=PVPX_TMP+V_TMP*(Y11-Y33)
     PVPY_TMP=PVPY_TMP+V_TMP*(X33-X11)
   END DO

   PUPX_TMP=PUPX_TMP/ART1(I)
   PUPY_TMP=PUPY_TMP/ART1(I)
   PVPX_TMP=PVPX_TMP/ART1(I)
   PVPY_TMP=PVPY_TMP/ART1(I)
   TMP1=PUPX_TMP**2+PVPY_TMP**2
   TMP2=0.5_SP*(PUPY_TMP+PVPX_TMP)**2
   VISCOFF(I)=SQRT(TMP1+TMP2)*ART1(I)

   IF(K == KBM1) THEN
     AH_BOTTOM(I) = (FACT*VISCOFF(I) + FM1)*NN_HVC(I)
   END IF

   DO II=1,NPCV
     I = NCEDGE_LST(II)
     I1=NTRG(I)
     IA=NIEC(I,1)
     IB=NIEC(I,2)
     
     IF((IA <= M .AND. IB <= M) .AND. I1 <= N)THEN
!       XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
!       YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
!!  ggao edge calculation
       XIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
       YIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))

       XIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       YIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       XI_TMP =0.5_SP*(XIJE1_TMP+XIJE2_TMP)
       YI_TMP =0.5_SP*(YIJE1_TMP+YIJE2_TMP)


       IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
!         XI_TMP = REARTH * COS(YI*DEG2RAD) * COS(XI*DEG2RAD) &
!                  * 2._SP /(1._SP+sin(YI*DEG2RAD))
!         YI_TMP = REARTH * COS(YI*DEG2RAD) * SIN(XI*DEG2RAD) &
!                  * 2._SP /(1._SP+sin(YI*DEG2RAD))

         VXA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * COS(VX(IA)*DEG2RAD) &
                   * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))
         VYA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * SIN(VX(IA)*DEG2RAD) &
                   * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))

         VXB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * COS(VX(IB)*DEG2RAD) &
                   * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))
         VYB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * SIN(VX(IB)*DEG2RAD) &
                   * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))

!         IF(IA == NODE_NORTHPOLE)THEN
         DXA=XI_TMP-VXA_TMP
         DYA=YI_TMP-VYA_TMP
!         ELSE IF(IB == NODE_NORTHPOLE)THEN
         DXB=XI_TMP-VXB_TMP
         DYB=YI_TMP-VYB_TMP
!	 END IF
!       END IF

        IF(IA == NODE_NORTHPOLE)THEN
	  PSPX_TMP=-PSPY(IB)*COS(VX(IB)*DEG2RAD)-PSPX(IB)*SIN(VX(IB)*DEG2RAD)
          PSPY_TMP=-PSPY(IB)*SIN(VX(IB)*DEG2RAD)+PSPX(IB)*COS(VX(IB)*DEG2RAD)
   
	  PSPXD_TMP=-PSPYD(IB)*COS(VX(IB)*DEG2RAD)-PSPXD(IB)*SIN(VX(IB)*DEG2RAD)
          PSPYD_TMP=-PSPYD(IB)*SIN(VX(IB)*DEG2RAD)+PSPXD(IB)*COS(VX(IB)*DEG2RAD)
   
          FIJ1=S1(IA,K)+DXA*PSPX(IA)+DYA*PSPY(IA)
          FIJ2=S1(IB,K)+DXB*PSPX_TMP+DYB*PSPY_TMP

          !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
          ! David moved HPRNU and added VHC
          VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) + FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB))) !/HPRNU

          TXX=0.5_SP*(PSPXD(IA)+PSPXD_TMP)*VISCOF
          TYY=0.5_SP*(PSPYD(IA)+PSPYD_TMP)*VISCOF
        ELSE IF(IB == NODE_NORTHPOLE)THEN
	  PSPX_TMP=-PSPY(IA)*COS(VX(IA)*DEG2RAD)-PSPX(IA)*SIN(VX(IA)*DEG2RAD)
          PSPY_TMP=-PSPY(IA)*SIN(VX(IA)*DEG2RAD)+PSPX(IA)*COS(VX(IA)*DEG2RAD)
   
	  PSPXD_TMP=-PSPYD(IA)*COS(VX(IA)*DEG2RAD)-PSPXD(IA)*SIN(VX(IA)*DEG2RAD)
          PSPYD_TMP=-PSPYD(IA)*SIN(VX(IA)*DEG2RAD)+PSPXD(IA)*COS(VX(IA)*DEG2RAD)
   
          FIJ1=S1(IA,K)+DXA*PSPX_TMP+DYA*PSPY_TMP
          FIJ2=S1(IB,K)+DXB*PSPX(IB)+DYB*PSPY(IB)

          !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
          ! David moved HPRNU and added VHC
          VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) + FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB)))

          TXX=0.5_SP*(PSPXD_TMP+PSPXD(IB))*VISCOF
          TYY=0.5_SP*(PSPYD_TMP+PSPYD(IB))*VISCOF
        END IF

       S1MIN=MINVAL(S1(NBSN(IA,1:NTSN(IA)-1),K))
       S1MIN=MIN(S1MIN, S1(IA,K))
       S1MAX=MAXVAL(S1(NBSN(IA,1:NTSN(IA)-1),K))
       S1MAX=MAX(S1MAX, S1(IA,K))
       S2MIN=MINVAL(S1(NBSN(IB,1:NTSN(IB)-1),K))
       S2MIN=MIN(S2MIN, S1(IB,K))
       S2MAX=MAXVAL(S1(NBSN(IB,1:NTSN(IB)-1),K))
       S2MAX=MAX(S2MAX, S1(IB,K))
       IF(FIJ1 < S1MIN) FIJ1=S1MIN
       IF(FIJ1 > S1MAX) FIJ1=S1MAX
       IF(FIJ2 < S2MIN) FIJ2=S2MIN
       IF(FIJ2 > S2MAX) FIJ2=S2MAX

!       IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
         UIJ_TMP = -V(I1,K)*COS(XC(I1)*DEG2RAD)-U(I1,K)*SIN(XC(I1)*DEG2RAD)
         VIJ_TMP = -V(I1,K)*SIN(XC(I1)*DEG2RAD)+U(I1,K)*COS(XC(I1)*DEG2RAD)
#        if defined (SEMI_IMPLICIT)
         UIJ_TMP1 = -VF(I1,K)*COS(XC(I1)*DEG2RAD)-UF(I1,K)*SIN(XC(I1)*DEG2RAD)
         VIJ_TMP1 = -VF(I1,K)*SIN(XC(I1)*DEG2RAD)+UF(I1,K)*COS(XC(I1)*DEG2RAD)
#        endif
       
         VX1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD)
         VY1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD)

         VX2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD)
         VY2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD)

         DLTXE_TMP = VX2_TMP-VX1_TMP
         DLTYE_TMP = VY2_TMP-VY1_TMP
       
         FXX=-DTIJ(I,K)*TXX*DLTYE_TMP
         FYY= DTIJ(I,K)*TYY*DLTXE_TMP

#        if !defined (SEMI_IMPLICIT)
         UVN_TMP = VIJ_TMP*DLTXE_TMP - UIJ_TMP*DLTYE_TMP  
         EXFLUX_TMP = -UVN_TMP*DTIJ(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP))*FIJ1)*0.5_SP
#        else
         UVN_TMP = VIJ_TMP*DLTXE_TMP - UIJ_TMP*DLTYE_TMP
         UVN_TMP1 = VIJ_TMP1*DLTXE_TMP - UIJ_TMP1*DLTYE_TMP
         EXFLUX_TMP = -UVN_TMP*DTIJ(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP))*FIJ1)*0.5_SP
         EXFLUX_TMP = (1.0_SP-CETA)*EXFLUX_TMP+CETA*( -UVN_TMP1*DTIJ1(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP1))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP1))*FIJ1)*0.5_SP )
#        endif       

         IF(IA == NODE_NORTHPOLE)THEN
           XFLUX(IA,K)=XFLUX(IA,K)+EXFLUX_TMP+FXX+FYY
           XFLUX_ADV(IA,K)=XFLUX_ADV(IA,K)+EXFLUX_TMP
         ELSE IF(IB == NODE_NORTHPOLE)THEN
           XFLUX(IB,K)=XFLUX(IB,K)-EXFLUX_TMP-FXX-FYY
           XFLUX_ADV(IB,K)=XFLUX_ADV(IB,K)-EXFLUX_TMP
         END IF
       END IF
     END IF  
   END DO

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: ADV_S_XY(K):",K

   RETURN
   END SUBROUTINE ADV_S_XY
!==============================================================================|

!==============================================================================|
   SUBROUTINE ADV_T_XY(XFLUX,XFLUX_ADV,PTPX,PTPY,PTPXD,PTPYD,VISCOFF,K,CETA)               

!------------------------------------------------------------------------------|

   IMPLICIT NONE
   INTEGER, INTENT(IN) :: K
   REAL(SP), DIMENSION(0:MT,KB)     :: XFLUX,XFLUX_ADV
   REAL(SP), DIMENSION(M)           :: PTPX,PTPY,PTPXD,PTPYD,VISCOFF
   REAL(SP), DIMENSION(3*(NT),KBM1)      :: DTIJ 
   REAL(SP) :: XI,YI
   REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2
   REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF
   REAL(SP) :: FACT,FM1
   INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,JTMP,JJ,II
   REAL(SP) :: TXPI,TYPI

   REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
   REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
   REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
   REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
   REAL(SP) :: PTPX_TMP,PTPY_TMP,PTPXD_TMP,PTPYD_TMP
   REAL(SP) :: U_TMP,V_TMP
   REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2

!!  ggao edge calculation
   REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
   REAL(SP) :: T1MIN, T1MAX, T2MIN, T2MAX

   REAL(SP) :: CETA
#  if defined (SEMI_IMPLICIT)
   REAL(SP), DIMENSION(3*(NT),KBM1)  :: DTIJ1
   REAL(SP) :: UIJ_TMP1, VIJ_TMP1, UVN_TMP1
#  endif
!------------------------------------------------------------------------------!

   ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN
   
   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: ADV_T_XY(K):",K

   
  SELECT CASE(HORIZONTAL_MIXING_TYPE)
  CASE ('closure')
     FACT = 1.0_SP
     FM1  = 0.0_SP
  CASE('constant')
     FACT = 0.0_SP
     FM1  = 1.0_SP
  CASE DEFAULT
     CALL FATAL_ERROR("UNKNOW HORIZONTAL MIXING TYPE:",&
          & TRIM(HORIZONTAL_MIXING_TYPE) )
  END SELECT

!
!--Initialize Fluxes-----------------------------------------------------------!
!
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     IA = NIEC(I,1)
     IB = NIEC(I,2)
     IF(IA == NODE_NORTHPOLE)THEN
       XFLUX(IA,K) = 0.0_SP
       XFLUX_ADV(IA,K) = 0.0_SP
     ELSE IF(IB == NODE_NORTHPOLE)THEN  
       XFLUX(IB,K) = 0.0_SP
       XFLUX_ADV(IB,K) = 0.0_SP
     END IF  
   END DO  
     
!
!--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
!
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     I1=NTRG(I)
     DTIJ(I,K)=DT1(I1)*DZ1(I1,K)
#    if defined (SEMI_IMPLICIT)
     DTIJ1(I,K) = D1(I1)*DZ1(I1,K)
#    endif
   END DO

!
!--Calculate the Advection and Horizontal Diffusion Terms----------------------!
!
   I = NODE_NORTHPOLE

   PUPX_TMP=0.0_SP
   PUPY_TMP=0.0_SP
   PVPX_TMP=0.0_SP
   PVPY_TMP=0.0_SP

   DO J=1,NTVE(I)
     I1=NBVE(I,J)
     JTMP=NBVT(I,J)
     J1=JTMP+1-(JTMP+1)/4*3
     J2=JTMP+2-(JTMP+2)/4*3
       
     VX_TMP = REARTH * COS(VY(I)*DEG2RAD) * COS(VX(I)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(I)*DEG2RAD))
     VY_TMP = REARTH * COS(VY(I)*DEG2RAD) * SIN(VX(I)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(I)*DEG2RAD))
		     
     VX1_TMP= REARTH * COS(VY(NV(I1,J1))*DEG2RAD) * COS(VX(NV(I1,J1))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J1))*DEG2RAD))
     VY1_TMP= REARTH * COS(VY(NV(I1,J1))*DEG2RAD) * SIN(VX(NV(I1,J1))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J1))*DEG2RAD))
		     
     VX2_TMP= REARTH * COS(YC(I1)*DEG2RAD) * COS(XC(I1)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(YC(I1)*DEG2RAD))
     VY2_TMP= REARTH * COS(YC(I1)*DEG2RAD) * SIN(XC(I1)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(YC(I1)*DEG2RAD))
		     
     VX3_TMP= REARTH * COS(VY(NV(I1,J2))*DEG2RAD) * COS(VX(NV(I1,J2))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J2))*DEG2RAD))
     VY3_TMP= REARTH * COS(VY(NV(I1,J2))*DEG2RAD) * SIN(VX(NV(I1,J2))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J2))*DEG2RAD))
		     
     X11=0.5_SP*(VX_TMP+VX1_TMP)
     Y11=0.5_SP*(VY_TMP+VY1_TMP)
     X22=VX2_TMP
     Y22=VX2_TMP
     X33=0.5_SP*(VX_TMP+VX3_TMP)
     Y33=0.5_SP*(VY_TMP+VY3_TMP)
     
     U_TMP = -V(I1,K)*COS(XC(I1)*DEG2RAD)-U(I1,K)*SIN(XC(I1)*DEG2RAD)
     V_TMP = -V(I1,K)*SIN(XC(I1)*DEG2RAD)+U(I1,K)*COS(XC(I1)*DEG2RAD)

     PUPX_TMP=PUPX_TMP+U_TMP*(Y11-Y33)
     PUPY_TMP=PUPY_TMP+U_TMP*(X33-X11)
     PVPX_TMP=PVPX_TMP+V_TMP*(Y11-Y33)
     PVPY_TMP=PVPY_TMP+V_TMP*(X33-X11)
   END DO

   PUPX_TMP=PUPX_TMP/ART1(I)
   PUPY_TMP=PUPY_TMP/ART1(I)
   PVPX_TMP=PVPX_TMP/ART1(I)
   PVPY_TMP=PVPY_TMP/ART1(I)
   TMP1=PUPX_TMP**2+PVPY_TMP**2
   TMP2=0.5_SP*(PUPY_TMP+PVPX_TMP)**2
   VISCOFF(I)=SQRT(TMP1+TMP2)*ART1(I)

   IF(K == KBM1) THEN
     AH_BOTTOM(I) = (FACT*VISCOFF(I) + FM1)*NN_HVC(I)
   END IF

   DO II=1,NPCV
     I = NCEDGE_LST(II)
     I1=NTRG(I)
     IA=NIEC(I,1)
     IB=NIEC(I,2)
     IF(IA <= M .AND. IB <= M .AND. I1 <= N)THEN
!       XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
!       YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
!!  ggao edge calculation
       XIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
       YIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))

       XIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       YIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       XI_TMP =0.5_SP*(XIJE1_TMP+XIJE2_TMP)
       YI_TMP =0.5_SP*(YIJE1_TMP+YIJE2_TMP)

       IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
!         XI_TMP = REARTH * COS(YI*DEG2RAD) * COS(XI*DEG2RAD) &
!                  * 2._SP /(1._SP+sin(YI*DEG2RAD))
!         YI_TMP = REARTH * COS(YI*DEG2RAD) * SIN(XI*DEG2RAD) &
!                  * 2._SP /(1._SP+sin(YI*DEG2RAD))

         VXA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * COS(VX(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))
         VYA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * SIN(VX(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))

         VXB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * COS(VX(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))
         VYB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * SIN(VX(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))

!        IF(IA == NODE_NORTHPOLE)THEN
           DXA=XI_TMP-VXA_TMP
           DYA=YI_TMP-VYA_TMP
!	 ELSE IF(IB == NODE_NORTHPOLE)THEN
           DXB=XI_TMP-VXB_TMP
           DYB=YI_TMP-VYB_TMP
!	 END IF
!       END IF

        IF(IA == NODE_NORTHPOLE)THEN
	  PTPX_TMP=-PTPY(IB)*COS(VX(IB)*DEG2RAD)-PTPX(IB)*SIN(VX(IB)*DEG2RAD)
          PTPY_TMP=-PTPY(IB)*SIN(VX(IB)*DEG2RAD)+PTPX(IB)*COS(VX(IB)*DEG2RAD)
   
	  PTPXD_TMP=-PTPYD(IB)*COS(VX(IB)*DEG2RAD)-PTPXD(IB)*SIN(VX(IB)*DEG2RAD)
          PTPYD_TMP=-PTPYD(IB)*SIN(VX(IB)*DEG2RAD)+PTPXD(IB)*COS(VX(IB)*DEG2RAD)
   
          FIJ1=T1(IA,K)+DXA*PTPX(IA)+DYA*PTPY(IA)
          FIJ2=T1(IB,K)+DXB*PTPX_TMP+DYB*PTPY_TMP

          !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
          ! David moved HPRNU and added VHC
          VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) + FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB))) !/HPRNU

          TXX=0.5_SP*(PTPXD(IA)+PTPXD_TMP)*VISCOF
          TYY=0.5_SP*(PTPYD(IA)+PTPYD_TMP)*VISCOF
        ELSE IF(IB == NODE_NORTHPOLE)THEN
	  PTPX_TMP=-PTPY(IA)*COS(VX(IA)*DEG2RAD)-PTPX(IA)*SIN(VX(IA)*DEG2RAD)
          PTPY_TMP=-PTPY(IA)*SIN(VX(IA)*DEG2RAD)+PTPX(IA)*COS(VX(IA)*DEG2RAD)
   
	  PTPXD_TMP=-PTPYD(IA)*COS(VX(IA)*DEG2RAD)-PTPXD(IA)*SIN(VX(IA)*DEG2RAD)
          PTPYD_TMP=-PTPYD(IA)*SIN(VX(IA)*DEG2RAD)+PTPXD(IA)*COS(VX(IA)*DEG2RAD)
   
          FIJ1=T1(IA,K)+DXA*PTPX_TMP+DYA*PTPY_TMP
          FIJ2=T1(IB,K)+DXB*PTPX(IB)+DYB*PTPY(IB)

          !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
          ! David moved HPRNU and added VHC
          VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) + FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB)))

          TXX=0.5_SP*(PTPXD_TMP+PTPXD(IB))*VISCOF
          TYY=0.5_SP*(PTPYD_TMP+PTPYD(IB))*VISCOF
        END IF

       T1MIN=MINVAL(T1(NBSN(IA,1:NTSN(IA)-1),K))
       T1MIN=MIN(T1MIN, T1(IA,K))
       T1MAX=MAXVAL(T1(NBSN(IA,1:NTSN(IA)-1),K))
       T1MAX=MAX(T1MAX, T1(IA,K))
       T2MIN=MINVAL(T1(NBSN(IB,1:NTSN(IB)-1),K))
       T2MIN=MIN(T2MIN, T1(IB,K))
       T2MAX=MAXVAL(T1(NBSN(IB,1:NTSN(IB)-1),K))
       T2MAX=MAX(T2MAX, T1(IB,K))
       IF(FIJ1 < T1MIN) FIJ1=T1MIN
       IF(FIJ1 > T1MAX) FIJ1=T1MAX
       IF(FIJ2 < T2MIN) FIJ2=T2MIN
       IF(FIJ2 > T2MAX) FIJ2=T2MAX

!     IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
         UIJ_TMP = -V(I1,K)*COS(XC(I1)*DEG2RAD)-U(I1,K)*SIN(XC(I1)*DEG2RAD)
         VIJ_TMP = -V(I1,K)*SIN(XC(I1)*DEG2RAD)+U(I1,K)*COS(XC(I1)*DEG2RAD)
#        if defined (SEMI_IMPLICIT)
         UIJ_TMP1 = -VF(I1,K)*COS(XC(I1)*DEG2RAD)-UF(I1,K)*SIN(XC(I1)*DEG2RAD)
         VIJ_TMP1 = -VF(I1,K)*SIN(XC(I1)*DEG2RAD)+UF(I1,K)*COS(XC(I1)*DEG2RAD)
#        endif
       
         VX1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD)
         VY1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD)

         VX2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD)
         VY2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD)

         DLTXE_TMP = VX2_TMP-VX1_TMP
         DLTYE_TMP = VY2_TMP-VY1_TMP
       
         FXX=-DTIJ(I,K)*TXX*DLTYE_TMP
         FYY= DTIJ(I,K)*TYY*DLTXE_TMP

#        if !defined (SEMI_IMPLICIT)
         UVN_TMP = VIJ_TMP*DLTXE_TMP - UIJ_TMP*DLTYE_TMP
         EXFLUX_TMP = -UVN_TMP*DTIJ(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP))*FIJ1)*0.5_SP
#        else
         UVN_TMP = VIJ_TMP*DLTXE_TMP - UIJ_TMP*DLTYE_TMP
         UVN_TMP1 = VIJ_TMP1*DLTXE_TMP - UIJ_TMP1*DLTYE_TMP
         EXFLUX_TMP = -UVN_TMP*DTIJ(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP))*FIJ1)*0.5_SP
         EXFLUX_TMP = (1.0_SP-CETA)*EXFLUX_TMP+CETA*( -UVN_TMP1*DTIJ1(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP1))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP1))*FIJ1)*0.5_SP )
#        endif
       
         IF(IA == NODE_NORTHPOLE)THEN
           XFLUX(IA,K)=XFLUX(IA,K)+EXFLUX_TMP+FXX+FYY
           XFLUX_ADV(IA,K)=XFLUX_ADV(IA,K)+EXFLUX_TMP
         ELSE IF(IB == NODE_NORTHPOLE)THEN
           XFLUX(IB,K)=XFLUX(IB,K)-EXFLUX_TMP-FXX-FYY
           XFLUX_ADV(IB,K)=XFLUX_ADV(IB,K)-EXFLUX_TMP
         END IF
       END IF
     END IF  
   END DO 

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: ADV_T_XY(K):",K

   RETURN
   END SUBROUTINE ADV_T_XY
!==============================================================================|

   SUBROUTINE ADV_N_XY(XFLUX,PWPX,PWPY,ISS,ID,DEP2,SPCDIR,N32)
!  This subroutine is for wave advection at the north pole    yzhang
!------------------------------------------------------------------------------|
   USE SWCOMM3, ONLY :MDC,MSC

   IMPLICIT NONE

   SAVE

   REAL :: N32(MDC,MSC,0:MT),N32_TMP(MDC,MSC,0:MT)
   INTEGER  :: ID,ISS,IG,IG2,IDT,IDD ! LWU IG2 IS FOR PLBC
   REAL :: CANX,CANY,CANX_TMP,CANY_TMP,ADDEXFLUX2455
   REAL :: SPCDIR(MDC,6)
   REAL(SP) :: DEP2(MT)
   REAL    :: DEPLOC,KWAVELOC,CGLOC,NN,ND,SPCSIGL
#  if defined (SPHERICAL)
   REAL(8) :: XTMP1,XTMP
#  endif
   REAL(SP), DIMENSION(0:MT)     :: XFLUX
   REAL(SP), DIMENSION(M)           :: PWPX,PWPY
   REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2
   INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,JTMP,JJ,II,L
   REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
   REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
   REAL(SP) :: UL_DEGREE,DL_DEGREE,CENTER_DEGREE,FF11
   REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
   REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
   REAL(SP) :: PWPX_TMP,PWPY_TMP
   REAL(SP) :: U_TMP,V_TMP
   REAL(SP) :: ADDYIJE1,ADDYIJE2
   REAL(SP) ::  ADD_DLTXE ,ADD_DLTYE,ADD_DLTXTRIE,ADD_DLTYTRIE
   REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2
   REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
   REAL(SP) :: AC1MIN, AC1MAX, AC2MIN, AC2MAX

#  if defined (SEMI_IMPLICIT)
   REAL(SP) :: UIJ_TMP1, VIJ_TMP1, UVN_TMP1
#  endif
!------------------------------------------------------------------------------!
   IF (NODE_NORTHPOLE .EQ. 0) RETURN

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: ADV_N_XY"

!   CAX = 0.0
!   CAY = 0.0

!
!--Initialize Fluxes-----------------------------------------------------------!
!
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     IA = NIEC(I,1)
     IB = NIEC(I,2)
     IF(IA == NODE_NORTHPOLE)THEN
       XFLUX(IA) = 0.0_SP
     ELSE IF(IB == NODE_NORTHPOLE)THEN
       XFLUX(IB) = 0.0_SP
     END IF
   END DO
!==========YZHANG_NORTHPOLE_TESTING================================
   I=NODE_NORTHPOLE
   PWPX(I)=0.0_SP
   PWPY(I)=0.0_SP
   ADD_DLTXTRIE=0.0_SP
   ADD_DLTYTRIE=0.0_SP
   DO J=1,NTSN(I)-1
     I1=NBSN(I,J)
     I2=NBSN(I,J+1)
     UL_DEGREE=22.5
     CENTER_DEGREE=0
     DL_DEGREE=337.5

     IF (VX(I1)<UL_DEGREE .OR. VX(I1)>DL_DEGREE)THEN
       IDT=MDC-(MDC*2/8)
       IDD=ID+IDT
!       IDD=ID
       IF(IDD>MDC)THEN
         IDD=IDD-MDC
       END IF
       N32_TMP(ID,ISS,I1)=N32(IDD,ISS,I1)
     END IF

     IF (VX(I2)<UL_DEGREE .OR. VX(I2)>DL_DEGREE)THEN
       IDT=MDC-(MDC*2/8)
       IDD=ID+IDT
!       IDD=ID
       IF(IDD>MDC)THEN
         IDD=IDD-MDC
       END IF
       N32_TMP(ID,ISS,I2)=N32(IDD,ISS,I2)
     END IF
     DO L=1,7
       CENTER_DEGREE=L*45.0
       UL_DEGREE=CENTER_DEGREE+22.5
       DL_DEGREE=CENTER_DEGREE-22.5

       IF(VX(I1)<UL_DEGREE .AND. VX(I1)>DL_DEGREE) THEN
         IDT=MDC-(MDC*(L+2)/8)
!         IDT=MDC-(MDC*L/8)
         IF(IDT<0)THEN
           IDT=IDT+MDC
         END IF
         IDD=ID+IDT
         IF(IDD>MDC)THEN
           IDD=IDD-MDC
         END IF
         N32_TMP(ID,ISS,I1)=N32(IDD,ISS,I1)
       END IF
       IF(VX(I2)<UL_DEGREE .AND. VX(I2)>DL_DEGREE) THEN
         IDT=MDC-(MDC*(L+2)/8)
!         IDT=MDC-(MDC*L/8)
         IF(IDT<0)THEN
           IDT=IDT+MDC
         END IF
         IDD=ID+IDT
         IF(IDD>MDC)THEN
           IDD=IDD-MDC
         END IF
         N32_TMP(ID,ISS,I2)=N32(IDD,ISS,I2)
       END IF
     END DO
     FF11=0.5*(N32_TMP(ID,ISS,I1)+N32_TMP(ID,ISS,I2))
     PWPX(I)=PWPX(I)+FF11*DLTYTRIE(I,J)
     PWPY(I)=PWPY(I)+FF11*DLTXTRIE(I,J)
     ADD_DLTXTRIE=ADD_DLTXTRIE+DLTXTRIE(I,J)
     ADD_DLTYTRIE=ADD_DLTYTRIE+DLTYTRIE(I,J)
   END DO

   PWPX(I)=PWPX(I)/ART2(I)
   PWPY(I)=PWPY(I)/ART2(I)
!============================================================================
   ADDYIJE1=0.0_SP
   ADDYIJE2=0.0_SP
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     ADDYIJE1=ADDYIJE1+YIJE(I,1)
     ADDYIJE2=ADDYIJE2+YIJE(I,2)
   END DO
   ADDYIJE1=ADDYIJE1/NPCV
   ADDYIJE2=ADDYIJE2/NPCV

   ADD_DLTXE=0.0_SP
   ADD_DLTYE=0.0_SP

   DO II=1,NPCV
     I = NCEDGE_LST(II)
     I1=NTRG(I)
     IA=NIEC(I,1)
     IB=NIEC(I,2)
     IF(IA <= M .AND. IB <= M .AND. I1 <= N)THEN
!        XI=0.5_SP*(XIJE(I,1)+XIJE(I,2))
!        YI=0.5_SP*(YIJE(I,1)+YIJE(I,2))
!!   ggao edge calculation
       XIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
       YIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))

       XIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       YIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       XI_TMP =0.5_SP*(XIJE1_TMP+XIJE2_TMP)
       YI_TMP =0.5_SP*(YIJE1_TMP+YIJE2_TMP)

       IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
         VXA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * COS(VX(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))
         VYA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * SIN(VX(IA)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))

         VXB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * COS(VX(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))
         VYB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * SIN(VX(IB)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))

!         IF(IA == NODE_NORTHPOLE)THEN
         DXA=XI_TMP-VXA_TMP
         DYA=YI_TMP-VYA_TMP
!         ELSE IF(IB == NODE_NORTHPOLE)THEN
         DXB=XI_TMP-VXB_TMP
         DYB=YI_TMP-VYB_TMP

         IF(IA == NODE_NORTHPOLE)THEN
           PWPX_TMP=-PWPY(IB)*COS(VX(IB)*DEG2RAD)-PWPX(IB)*SIN(VX(IB)*DEG2RAD)
           PWPY_TMP=-PWPY(IB)*SIN(VX(IB)*DEG2RAD)+PWPX(IB)*COS(VX(IB)*DEG2RAD)

!           PSPXD_TMP=-PSPYD(IB)*COS(VX(IB)*DEG2RAD)-PSPXD(IB)*SIN(VX(IB)*DEG2RAD)
!           PSPYD_TMP=-PSPYD(IB)*SIN(VX(IB)*DEG2RAD)+PSPXD(IB)*COS(VX(IB)*DEG2RAD)
!======================zhangyang======start=======================
           IF (VX(IB)<UL_DEGREE .OR. VX(IB)>DL_DEGREE)THEN
             IDT=MDC-(MDC*2/8)
             IDD=ID+IDT
!             IDD=ID
             IF(IDD>MDC)THEN
               IDD=IDD-MDC
             END IF
             N32_TMP(ID,ISS,IB)=N32(IDD,ISS,IB)
           END IF

           DO L=1,7
             CENTER_DEGREE=L*45.0
             UL_DEGREE=CENTER_DEGREE+22.5
             DL_DEGREE=CENTER_DEGREE-22.5
             IF(VX(IB)<UL_DEGREE .AND. VX(IB)>DL_DEGREE) THEN
               IDT=MDC-(MDC*(L+2)/8)
!               IDT=MDC-(MDC*L/8)
               IF(IDT<0)THEN
                 IDT=IDT+MDC
               END IF
               IDD=ID+IDT
               IF(IDD>MDC)THEN
                 IDD=IDD-MDC
               END IF
               N32_TMP(ID,ISS,IB)=N32(IDD,ISS,IB)
             END IF
           END DO
           FIJ1=N32(ID,ISS,IA)!+DXA*PWPX(IA)+DYA*PWPY(IA)
           FIJ2=N32_TMP(ID,ISS,IB)!+DXB*PWPX_TMP+DYB*PWPY_TMP
!================zhangyang======end==================================

          !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
          ! David moved HPRNU and added VHC
!          VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) +
!          FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB))) !/HPRNU

!          TXX=0.5_SP*(PSPXD(IA)+PSPXD_TMP)*VISCOF
!          TYY=0.5_SP*(PSPYD(IA)+PSPYD_TMP)*VISCOF

         ELSE IF(IB == NODE_NORTHPOLE)THEN
           PWPX_TMP=-PWPY(IA)*COS(VX(IA)*DEG2RAD)-PWPX(IA)*SIN(VX(IA)*DEG2RAD)
           PWPY_TMP=-PWPY(IA)*SIN(VX(IA)*DEG2RAD)+PWPX(IA)*COS(VX(IA)*DEG2RAD)

!           PSPXD_TMP=-PSPYD(IA)*COS(VX(IA)*DEG2RAD)-PSPXD(IA)*SIN(VX(IA)*DEG2RAD)
!           PSPYD_TMP=-PSPYD(IA)*SIN(VX(IA)*DEG2RAD)+PSPXD(IA)*COS(VX(IA)*DEG2RAD)

!======================zhangyang======start=======================
           IF (VX(IA)<UL_DEGREE .OR. VX(IA)>DL_DEGREE)THEN
             IDT=MDC-(MDC*2/8)
             IDD=ID+IDT
!             IDD=ID
             IF(IDD>MDC)THEN
               IDD=IDD-MDC
             END IF
             N32_TMP(ID,ISS,IA)=N32(IDD,ISS,IA)
           END IF

           DO L=1,7
             CENTER_DEGREE=L*45.0
             UL_DEGREE=CENTER_DEGREE+22.5
             DL_DEGREE=CENTER_DEGREE-22.5
             IF(VX(IA)<UL_DEGREE .AND. VX(IA)>DL_DEGREE) THEN
               IDT=MDC-(MDC*(L+2)/8)
!               IDT=MDC-(MDC*L/8)
               IF(IDT<0)THEN
                 IDT=IDT+MDC
               END IF
               IDD=ID+IDT
               IF(IDD>MDC)THEN
                 IDD=IDD-MDC
               END IF
               N32_TMP(ID,ISS,IA)=N32(IDD,ISS,IA)
             END IF
           END DO
           FIJ1=N32_TMP(ID,ISS,IA)!+DXA*PWPX_TMP+DYA*PWPY_TMP
           FIJ2=N32(ID,ISS,IB)!+DXB*PWPX(IB)+DYB*PWPY(IB)

!================zhangyang======end==================================
         END IF
         CALL SWAPAR1(I1,ISS,ID,DEP2(1),KWAVELOC,CGLOC)
         UIJ_TMP = UA(I1)
         VIJ_TMP = VA(I1)

         DO L=1,8
           CENTER_DEGREE=L*45.0-22.5
           UL_DEGREE=CENTER_DEGREE+22.5
           DL_DEGREE=CENTER_DEGREE-22.5

           IF(XC(I1)<UL_DEGREE .AND. XC(I1)>DL_DEGREE) THEN
             IDT=MDC-(MDC*(L+2)/8-5)
!             IDT=MDC-(MDC*L/8)
             IF(IDT<0)THEN
               IDT=IDT+MDC
             END IF
             IDD=ID+IDT
             IF(IDD>MDC)THEN
               IDD=IDD-MDC
             END IF
           END IF
         END DO

         CALL SPROXY(I1,ISS,IDD,CANX,CANY,CGLOC,SPCDIR(IDD,2),SPCDIR(IDD,3),UIJ_TMP,VIJ_TMP)

         VX1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
         VY1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
 
         VX2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
         VY2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))

         DLTXE_TMP = VX2_TMP-VX1_TMP
         DLTYE_TMP = VY2_TMP-VY1_TMP

         CANX_TMP = -CANY*COS(XC(I1)*DEG2RAD)-CANX*SIN(XC(I1)*DEG2RAD)
         CANY_TMP = -CANY*SIN(XC(I1)*DEG2RAD)+CANX*COS(XC(I1)*DEG2RAD)

#        if !defined (SEMI_IMPLICIT)
         UVN_TMP = CANY_TMP*DLTXE_TMP - CANX_TMP*DLTYE_TMP
         EXFLUX_TMP = -UVN_TMP*((1.0_SP+SIGN(1.0_SP,UVN_TMP))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP))*FIJ1)*0.5_SP
#        endif

         IF(IA == NODE_NORTHPOLE)THEN
           XFLUX(IA)=XFLUX(IA)+EXFLUX_TMP
!           XFLUX_ADV(IA,K)=XFLUX_ADV(IA,K)+EXFLUX_TMP
         ELSE IF(IB == NODE_NORTHPOLE)THEN
           XFLUX(IB)=XFLUX(IB)-EXFLUX_TMP
!           XFLUX_ADV(IB,K)=XFLUX_ADV(IB,K)-EXFLUX_TMP
         END IF
       END IF
     END IF
   END DO

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: ADV_N_XY(ID,ISS):",ID,ISS

   RETURN
   END SUBROUTINE ADV_N_XY

!==============================================================================|
!     CALCULATE THE BAROCLINIC PRESSURE GRADIENT IN SIGMA COORDINATES          |
!==============================================================================|

   SUBROUTINE BAROPG_XY(DRIJK1,DRIJK2) 

!==============================================================================|
   IMPLICIT NONE
   REAL(SP) :: DRIJK1(0:N,3,KBM1), DRIJK2(0:N,KBM1)
   REAL(SP) :: DIJ,DRHO1,DRHO2
   INTEGER  :: I,II,K,J,J1,J2,IJK
   REAL(SP) :: VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP
!==============================================================================|

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: baropg_XY"


   DO II = 1, NP
     I=NP_LST(II)
     DO K=1,KBM1
       DRHOX(I,K)=0.0_SP
       DRHOY(I,K)=0.0_SP
     END DO
   END DO

   DO II = 1, NP
     I=NP_LST(II)
     DO K=1,KBM1
        DO J = 1, 3
          J1=J+1-INT((J+1)/4)*3
          J2=J+2-INT((J+2)/4)*3
          IJK=NBE(I,J)
          DIJ=0.5_SP*(DT(NV(I,J1))+DT(NV(I,J2)))

          VY1_TMP=REARTH*COS(VY(NV(I,J1))*DEG2RAD)*SIN(VX(NV(I,J1))*DEG2RAD)
          VY2_TMP=REARTH*COS(VY(NV(I,J2))*DEG2RAD)*SIN(VX(NV(I,J2))*DEG2RAD)

          DRHO1=(VY1_TMP-VY2_TMP)*DRIJK1(I,J,K)*DT1(I)
          DRHO2=(VY1_TMP-VY2_TMP)*DIJ*DRIJK2(I,K)

          DRHOX(I,K)=DRHOX(I,K)+DRHO1+DRHO2

          VX1_TMP=REARTH*COS(VY(NV(I,J1))*DEG2RAD)*COS(VX(NV(I,J1))*DEG2RAD)
          VX2_TMP=REARTH*COS(VY(NV(I,J2))*DEG2RAD)*COS(VX(NV(I,J2))*DEG2RAD)

	  DRHO1=(VX2_TMP-VX1_TMP)*DRIJK1(I,J,K)*DT1(I)
          DRHO2=(VX2_TMP-VX1_TMP)*DIJ*DRIJK2(I,K)

          DRHOY(I,K)=DRHOY(I,K)+DRHO1+DRHO2

       END DO
     END DO
   END DO

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: baropg_XY"
   
   RETURN
   END SUBROUTINE BAROPG_XY
!==============================================================================|

!==============================================================================|
   SUBROUTINE SHAPE_COEF_XY

!----------------------------------------------------------------------!
!  This subrountine is used to calculate the coefficient for a linear  !
!  function on the x-y plane, i.e.:                                    !
!                     r(x,y;phai)=phai_c+cofa1*x+cofa2*y               !
!     innc(i)=0    cells on the boundary                               !
!     innc(i)=1    cells in the interior                               !
!----------------------------------------------------------------------!
     
   USE ALL_VARS
   IMPLICIT NONE
   REAL(DP) X1,X2,X3,Y1,Y2,Y3,DELT,AI1,AI2,AI3,BI1,BI2,BI3,CI1,CI2,CI3
   REAL(DP) DELTX,DELTY,TEMP1,ANG1,ANG2,B1,B2,ANGLE
   REAL(DP), ALLOCATABLE :: XC_TMP(:),YC_TMP(:),VX_TMP(:),VY_TMP(:)
   INTEGER  I,II,J,JJ,J1,J2
!
!---------------interior cells-----------------------------------------!
!
   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: shape_coef_xy"

   ! ALL OTHER PROCESSORS ESCAPE HERE
   IF (NODE_NORTHPOLE .EQ. 0) RETURN


   ALLOCATE(A1U_XY(N,4)); A1U_XY = 0.0_SP
   ALLOCATE(A2U_XY(N,4)); A2U_XY = 0.0_SP
   ALLOCATE(AW0_XY(N,3)); AW0_XY = 0.0_SP
   ALLOCATE(AWX_XY(N,3)); AWX_XY = 0.0_SP
   ALLOCATE(AWY_XY(N,3)); AWY_XY = 0.0_SP
   
   ALLOCATE(XC_TMP(0:NT)); XC_TMP = 0.0_SP
   ALLOCATE(YC_TMP(0:NT)); YC_TMP = 0.0_SP
   ALLOCATE(VX_TMP(0:MT)); VX_TMP = 0.0_SP
   ALLOCATE(VY_TMP(0:MT)); VY_TMP = 0.0_SP
   
   DO I=1,NT
     XC_TMP(I) = REARTH * COS(YC(I)*DEG2RAD) * COS(XC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YC(I)*DEG2RAD))
     YC_TMP(I) = REARTH * COS(YC(I)*DEG2RAD) * SIN(XC(I)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YC(I)*DEG2RAD))
   END DO		  

   DO I=1,MT
     VX_TMP(I) = REARTH * COS(VY(I)*DEG2RAD) * COS(VX(I)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(I)*DEG2RAD))
     VY_TMP(I) = REARTH * COS(VY(I)*DEG2RAD) * SIN(VX(I)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(VY(I)*DEG2RAD))
   END DO		  

   DO I=1,N
     IF(ISBCE(I) == 0)THEN
       Y1 = YC_TMP(NBE(I,1))-YC_TMP(I)
       Y2 = YC_TMP(NBE(I,2))-YC_TMP(I)
       Y3 = YC_TMP(NBE(I,3))-YC_TMP(I)
       X1=XC_TMP(NBE(I,1))-XC_TMP(I)
       X2=XC_TMP(NBE(I,2))-XC_TMP(I)
       X3=XC_TMP(NBE(I,3))-XC_TMP(I)

       X1=X1/1000.0_SP
       X2=X2/1000.0_SP
       X3=X3/1000.0_SP
       Y1=Y1/1000.0_SP
       Y2=Y2/1000.0_SP
       Y3=Y3/1000.0_SP

       delt=(x1*y2-x2*y1)**2+(x1*y3-x3*y1)**2+(x2*y3-x3*y2)**2
       delt=delt*1000.0_SP

       a1u_XY(i,1)=(y1+y2+y3)*(x1*y1+x2*y2+x3*y3)- &
                (x1+x2+x3)*(y1**2+y2**2+y3**2)
       a1u_XY(i,1)=a1u_XY(i,1)/delt
       a1u_XY(i,2)=(y1**2+y2**2+y3**2)*x1-(x1*y1+x2*y2+x3*y3)*y1
       a1u_XY(i,2)=a1u_XY(i,2)/delt
       a1u_XY(i,3)=(y1**2+y2**2+y3**2)*x2-(x1*y1+x2*y2+x3*y3)*y2
       a1u_XY(i,3)=a1u_XY(i,3)/delt
       a1u_XY(i,4)=(y1**2+y2**2+y3**2)*x3-(x1*y1+x2*y2+x3*y3)*y3
       a1u_XY(i,4)=a1u_XY(i,4)/delt

       a2u_XY(i,1)=(x1+x2+x3)*(x1*y1+x2*y2+x3*y3)- &
                (y1+y2+y3)*(x1**2+x2**2+x3**2)
       a2u_XY(i,1)=a2u_XY(i,1)/delt
       a2u_XY(i,2)=(x1**2+x2**2+x3**2)*y1-(x1*y1+x2*y2+x3*y3)*x1
       a2u_XY(i,2)=a2u_XY(i,2)/delt
       a2u_XY(i,3)=(x1**2+x2**2+x3**2)*y2-(x1*y1+x2*y2+x3*y3)*x2
       a2u_XY(i,3)=a2u_XY(i,3)/delt
       a2u_XY(i,4)=(x1**2+x2**2+x3**2)*y3-(x1*y1+x2*y2+x3*y3)*x3
       a2u_XY(i,4)=a2u_XY(i,4)/delt
     end if

     x1=vx_TMP(nv(i,1))-xc_TMP(i)
     x2=vx_TMP(nv(i,2))-xc_TMP(i)
     x3=vx_TMP(nv(i,3))-xc_TMP(i)
     y1=vy_TMP(nv(i,1))-yc_TMP(i)
     y2=vy_TMP(nv(i,2))-yc_TMP(i)
     y3=vy_TMP(nv(i,3))-yc_TMP(i)


     ai1=y2-y3
     ai2=y3-y1
     ai3=y1-y2
     bi1=x3-x2
     bi2=x1-x3
     bi3=x2-x1
     ci1=x2*y3-x3*y2
     ci2=x3*y1-x1*y3
     ci3=x1*y2-x2*y1

     aw0_XY(i,1)=-ci1/2./art(i)
     aw0_XY(i,2)=-ci2/2./art(i)
     aw0_XY(i,3)=-ci3/2./art(i)
     awx_XY(i,1)=-ai1/2./art(i)
     awx_XY(i,2)=-ai2/2./art(i)
     awx_XY(i,3)=-ai3/2./art(i)
     awy_XY(i,1)=-bi1/2./art(i)
     awy_XY(i,2)=-bi2/2./art(i)
     awy_XY(i,3)=-bi3/2./art(i)
   end do

   DEALLOCATE(XC_TMP,YC_TMP,VX_TMP,VY_TMP)
   
   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: shape_coef_xy"

   return
   end subroutine shape_coef_xy

!==============================================================================|

!==============================================================================|
!==============================================================================|
   SUBROUTINE ADV_Q_XY(XFLUX,PQPX,PQPY,PQPXD,PQPYD,VISCOFF,Q,UQ,VQ,K,UQ1,VQ1,CETA)               

!==============================================================================|
!   Calculate the Turbulent Kinetic Energy and Mixing Length Based on          |
!   The Mellor-Yamada Level 2.5 Turbulent Closure Model                        |
!==============================================================================|


!------------------------------------------------------------------------------|

   IMPLICIT NONE
   INTEGER, INTENT(IN) :: K
   REAL(SP), DIMENSION(0:MT,KB)     :: XFLUX,Q
   REAL(SP), DIMENSION(M)           :: PQPX,PQPY,PQPXD,PQPYD,VISCOFF
   REAL(SP), DIMENSION(3*(NT),KBM1)      :: DTIJ 
   REAL(SP) :: XI,YI
   REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2 
   REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF   
   REAL(SP) :: FACT,FM1
   INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,JTMP,JJ,II
   REAL(SP) :: TXPI,TYPI

   REAL(SP) :: VX_TMP,VY_TMP,VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP,VX3_TMP,VY3_TMP
   REAL(SP) :: XI_TMP,YI_TMP,VXA_TMP,VYA_TMP,VXB_TMP,VYB_TMP
   REAL(SP) :: UIJ_TMP,VIJ_TMP,DLTXE_TMP,DLTYE_TMP,UVN_TMP,EXFLUX_TMP
   REAL(SP) :: PUPX_TMP,PUPY_TMP,PVPX_TMP,PVPY_TMP
   REAL(SP) :: PQPX_TMP,PQPY_TMP,PQPXD_TMP,PQPYD_TMP
   REAL(SP) :: U_TMP,V_TMP
   REAL(SP) :: X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2

   REAL(SP) :: XIJE1_TMP,YIJE1_TMP,XIJE2_TMP,YIJE2_TMP
   REAL(SP) :: Q1MIN, Q1MAX, Q2MIN, Q2MAX

   REAL(SP), DIMENSION(0:NT,KB)    :: UQ,VQ

   REAL(SP), DIMENSION(0:,:)    :: UQ1, VQ1
   REAL(SP) :: CETA
#  if defined (SEMI_IMPLICIT)
   REAL(SP), DIMENSION(3*(NT),KBM1)  :: DTIJ1
   REAL(SP) :: UIJ_TMP1, VIJ_TMP1, UVN_TMP1
#  endif
!------------------------------------------------------------------------------!
   IF (NODE_NORTHPOLE .EQ. 0) RETURN

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: adv_q_xy(K):",K


   SELECT CASE(HORIZONTAL_MIXING_TYPE)
   CASE ('closure')
      FACT = 1.0_SP
      FM1  = 0.0_SP
   CASE('constant')
      FACT = 0.0_SP
      FM1  = 1.0_SP
   CASE DEFAULT
      CALL FATAL_ERROR("UNKNOW HORIZONTAL MIXING TYPE:",&
           & TRIM(HORIZONTAL_MIXING_TYPE) )
   END SELECT
   

!
!--Initialize Fluxes-----------------------------------------------------------!
!
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     IA = NIEC(I,1)
     IB = NIEC(I,2)
     IF(IA == NODE_NORTHPOLE)THEN
       XFLUX(IA,K) = 0.0_SP
     ELSE IF(IB == NODE_NORTHPOLE)THEN  
       XFLUX(IB,K) = 0.0_SP
     END IF  
   END DO  
     
!
!--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------!
!
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     I1=NTRG(I)
     DTIJ(I,K)=DT1(I1)*DZ1(I1,K)
#    if defined (SEMI_IMPLICIT)
     DTIJ1(I,K)=D1(I1)*DZ1(I1,K)
#    endif
   END DO

!
!--Calculate the Advection and Horizontal Diffusion Terms----------------------!
!
   I = NODE_NORTHPOLE

   IF(I==0)  RETURN

   IF(I /= 0)THEN
   PUPX_TMP=0.0_SP
   PUPY_TMP=0.0_SP
   PVPX_TMP=0.0_SP
   PVPY_TMP=0.0_SP

   DO J=1,NTVE(I)
     I1=NBVE(I,J)
     JTMP=NBVT(I,J)
     J1=JTMP+1-(JTMP+1)/4*3
     J2=JTMP+2-(JTMP+2)/4*3
      
     VX_TMP = REARTH * COS(VY(I)*DEG2RAD) * COS(VX(I)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(I)*DEG2RAD))
     VY_TMP = REARTH * COS(VY(I)*DEG2RAD) * SIN(VX(I)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(I)*DEG2RAD))
		     
     VX1_TMP= REARTH * COS(VY(NV(I1,J1))*DEG2RAD) * COS(VX(NV(I1,J1))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J1))*DEG2RAD))
     VY1_TMP= REARTH * COS(VY(NV(I1,J1))*DEG2RAD) * SIN(VX(NV(I1,J1))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J1))*DEG2RAD))
		     
     VX2_TMP= REARTH * COS(YC(I1)*DEG2RAD) * COS(XC(I1)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(YC(I1)*DEG2RAD))
     VY2_TMP= REARTH * COS(YC(I1)*DEG2RAD) * SIN(XC(I1)*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(YC(I1)*DEG2RAD))
		     
     VX3_TMP= REARTH * COS(VY(NV(I1,J2))*DEG2RAD) * COS(VX(NV(I1,J2))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J2))*DEG2RAD))
     VY3_TMP= REARTH * COS(VY(NV(I1,J2))*DEG2RAD) * SIN(VX(NV(I1,J2))*DEG2RAD) &
                     * 2._SP /(1._SP+SIN(VY(NV(I1,J2))*DEG2RAD))
		     
     X11=0.5_SP*(VX_TMP+VX1_TMP)
     Y11=0.5_SP*(VY_TMP+VY1_TMP)
     X22=VX2_TMP
     Y22=VX2_TMP
     X33=0.5_SP*(VX_TMP+VX3_TMP)
     Y33=0.5_SP*(VY_TMP+VY3_TMP)
     
     U_TMP = -VQ(I1,K)*COS(XC(I1)*DEG2RAD)-UQ(I1,K)*SIN(XC(I1)*DEG2RAD)
     V_TMP = -VQ(I1,K)*SIN(XC(I1)*DEG2RAD)+UQ(I1,K)*COS(XC(I1)*DEG2RAD)

     PUPX_TMP=PUPX_TMP+U_TMP*(Y11-Y33)
     PUPY_TMP=PUPY_TMP+U_TMP*(X33-X11)
     PVPX_TMP=PVPX_TMP+V_TMP*(Y11-Y33)
     PVPY_TMP=PVPY_TMP+V_TMP*(X33-X11)
   END DO

   PUPX_TMP=PUPX_TMP/ART1(I)
   PUPY_TMP=PUPY_TMP/ART1(I)
   PVPX_TMP=PVPX_TMP/ART1(I)
   PVPY_TMP=PVPY_TMP/ART1(I)
   TMP1=PUPX_TMP**2+PVPY_TMP**2
   TMP2=0.5_SP*(PUPY_TMP+PVPX_TMP)**2
   VISCOFF(I)=SQRT(TMP1+TMP2)*ART1(I)

   IF(K == KBM1) THEN
     AH_BOTTOM(I) = (FACT*VISCOFF(I) + FM1)*NN_HVC(I)
   END IF
   endif
   DO II=1,NPCV
     I = NCEDGE_LST(II)
     I1=NTRG(I)
     IA=NIEC(I,1)
     IB=NIEC(I,2)
     
     IF((IA <= M .AND. IB <= M) .AND. I1 <= N)THEN
       XIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))
       YIJE1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,1)*DEG2RAD))

       XIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       YIJE2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD) &
                  * 2._SP /(1._SP+sin(YIJE(I,2)*DEG2RAD))
       XI_TMP =0.5_SP*(XIJE1_TMP+XIJE2_TMP)
       YI_TMP =0.5_SP*(YIJE1_TMP+YIJE2_TMP)

       IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
         VXA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * COS(VX(IA)*DEG2RAD) &
                   * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))
         VYA_TMP = REARTH * COS(VY(IA)*DEG2RAD) * SIN(VX(IA)*DEG2RAD) &
                   * 2._SP /(1._SP+sin(VY(IA)*DEG2RAD))

         VXB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * COS(VX(IB)*DEG2RAD) &
                   * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))
         VYB_TMP = REARTH * COS(VY(IB)*DEG2RAD) * SIN(VX(IB)*DEG2RAD) &
                   * 2._SP /(1._SP+sin(VY(IB)*DEG2RAD))

         DXA=XI_TMP-VXA_TMP
         DYA=YI_TMP-VYA_TMP
         DXB=XI_TMP-VXB_TMP
         DYB=YI_TMP-VYB_TMP

        IF(IA == NODE_NORTHPOLE)THEN
	  PQPX_TMP=-PQPY(IB)*COS(VX(IB)*DEG2RAD)-PQPX(IB)*SIN(VX(IB)*DEG2RAD)
          PQPY_TMP=-PQPY(IB)*SIN(VX(IB)*DEG2RAD)+PQPX(IB)*COS(VX(IB)*DEG2RAD)
   
	  PQPXD_TMP=-PQPYD(IB)*COS(VX(IB)*DEG2RAD)-PQPXD(IB)*SIN(VX(IB)*DEG2RAD)
          PQPYD_TMP=-PQPYD(IB)*SIN(VX(IB)*DEG2RAD)+PQPXD(IB)*COS(VX(IB)*DEG2RAD)
   
          FIJ1=Q(IA,K)+DXA*PQPX(IA)+DYA*PQPY(IA)
          FIJ2=Q(IB,K)+DXB*PQPX_TMP+DYB*PQPY_TMP

          ! VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
          ! David moved HPRNU and added VHC
          VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) + FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB)))/HPRNU

          TXX=0.5_SP*(PQPXD(IA)+PQPXD_TMP)*VISCOF
          TYY=0.5_SP*(PQPYD(IA)+PQPYD_TMP)*VISCOF
        ELSE IF(IB == NODE_NORTHPOLE)THEN
	  PQPX_TMP=-PQPY(IA)*COS(VX(IA)*DEG2RAD)-PQPX(IA)*SIN(VX(IA)*DEG2RAD)
          PQPY_TMP=-PQPY(IA)*SIN(VX(IA)*DEG2RAD)+PQPX(IA)*COS(VX(IA)*DEG2RAD)
   
	  PQPXD_TMP=-PQPYD(IA)*COS(VX(IA)*DEG2RAD)-PQPXD(IA)*SIN(VX(IA)*DEG2RAD)
          PQPYD_TMP=-PQPYD(IA)*SIN(VX(IA)*DEG2RAD)+PQPXD(IA)*COS(VX(IA)*DEG2RAD)
   
          FIJ1=Q(IA,K)+DXA*PQPX_TMP+DYA*PQPY_TMP
          FIJ2=Q(IB,K)+DXB*PQPX(IB)+DYB*PQPY(IB)

          !VISCOF=HORCON*(FACT*(VISCOFF(IA)+VISCOFF(IB))*0.5_SP + FM1)
          ! David moved HPRNU and added VHC
          VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) + FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB)))/HPRNU

          TXX=0.5_SP*(PQPXD_TMP+PQPXD(IB))*VISCOF
          TYY=0.5_SP*(PQPYD_TMP+PQPYD(IB))*VISCOF
        END IF

       Q1MIN=MINVAL(Q(NBSN(IA,1:NTSN(IA)-1),K))
       Q1MIN=MIN(Q1MIN, Q(IA,K))
       Q1MAX=MAXVAL(Q(NBSN(IA,1:NTSN(IA)-1),K))
       Q1MAX=MAX(Q1MAX, Q(IA,K))
       Q2MIN=MINVAL(Q(NBSN(IB,1:NTSN(IB)-1),K))
       Q2MIN=MIN(Q2MIN, Q(IB,K))
       Q2MAX=MAXVAL(Q(NBSN(IB,1:NTSN(IB)-1),K))
       Q2MAX=MAX(Q2MAX, Q(IB,K))
       IF(FIJ1 < Q1MIN) FIJ1=Q1MIN
       IF(FIJ1 > Q1MAX) FIJ1=Q1MAX
       IF(FIJ2 < Q2MIN) FIJ2=Q2MIN
       IF(FIJ2 > Q2MAX) FIJ2=Q2MAX

!        FXX=-DTIJ(I,K)*TXX*DLTYE(I)
!        FYY= DTIJ(I,K)*TYY*DLTXE(I)

!       IF(IA == NODE_NORTHPOLE .OR. IB == NODE_NORTHPOLE)THEN
         UIJ_TMP = -VQ(I1,K)*COS(XC(I1)*DEG2RAD)-UQ(I1,K)*SIN(XC(I1)*DEG2RAD)
         VIJ_TMP = -VQ(I1,K)*SIN(XC(I1)*DEG2RAD)+UQ(I1,K)*COS(XC(I1)*DEG2RAD)
#        if defined (SEMI_IMPLICIT)
         UIJ_TMP1 = -VQ1(I1,K)*COS(XC(I1)*DEG2RAD)-UQ1(I1,K)*SIN(XC(I1)*DEG2RAD)
         VIJ_TMP1 = -VQ1(I1,K)*SIN(XC(I1)*DEG2RAD)+UQ1(I1,K)*COS(XC(I1)*DEG2RAD)
#        endif
       
         VX1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * COS(XIJE(I,1)*DEG2RAD)
         VY1_TMP = REARTH * COS(YIJE(I,1)*DEG2RAD) * SIN(XIJE(I,1)*DEG2RAD)

         VX2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * COS(XIJE(I,2)*DEG2RAD)
         VY2_TMP = REARTH * COS(YIJE(I,2)*DEG2RAD) * SIN(XIJE(I,2)*DEG2RAD)

         DLTXE_TMP = VX2_TMP-VX1_TMP
         DLTYE_TMP = VY2_TMP-VY1_TMP
       
         FXX=-DTIJ(I,K)*TXX*DLTYE_TMP
         FYY= DTIJ(I,K)*TYY*DLTXE_TMP

#        if !defined (SEMI_IMPLICIT)
         UVN_TMP = VIJ_TMP*DLTXE_TMP - UIJ_TMP*DLTYE_TMP
         EXFLUX_TMP = -UVN_TMP*DTIJ(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP))*FIJ1)*0.5_SP
#        else
         UVN_TMP = VIJ_TMP*DLTXE_TMP - UIJ_TMP*DLTYE_TMP
         UVN_TMP1 = VIJ_TMP1*DLTXE_TMP - UIJ_TMP1*DLTYE_TMP
         EXFLUX_TMP = -UVN_TMP*DTIJ(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP))*FIJ1)*0.5_SP
         EXFLUX_TMP = (1.0_SP-CETA)*EXFLUX_TMP+CETA*( -UVN_TMP1*DTIJ1(I,K)*((1.0_SP+SIGN(1.0_SP,UVN_TMP1))*FIJ2+   &
                      (1.0_SP-SIGN(1.0_SP,UVN_TMP1))*FIJ1)*0.5_SP )
#        endif
       
         IF(IA == NODE_NORTHPOLE)THEN
           XFLUX(IA,K)=XFLUX(IA,K)+EXFLUX_TMP+FXX+FYY
         ELSE IF(IB == NODE_NORTHPOLE)THEN
           XFLUX(IB,K)=XFLUX(IB,K)-EXFLUX_TMP-FXX-FYY
         END IF
       END IF
     END IF  
   END DO

   IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: ADV_Q_XY(K):",K

   RETURN
   END SUBROUTINE ADV_Q_XY
!==============================================================================|

END MODULE MOD_NORTHPOLE
